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O ; ABSTRACT 

• Radio galaxies and quasars are among the largest and most powerful single objects 

known and are believed to have had a significant impact on the evolving Universe and 

I its large scale structure. We explore the intrinsic and extrinsic properties of the popu- 

lation of Fanaroff-Riley II (FR II) objects, that is their kinetic luminosities, lifetimes, 
Q | and the central densities of their environments. In particular, the radio and kinetic lu- 

S—{ , minosity functions of these powerful radio sources are investigated using the complete, 

C/3 ■ flux limited radio catalogues of 3CRR and Best et al.. We construct multidimensional 

| Monte Carlo simulations using semi-analytical models of FR II source time evolution 

to create artificial samples of radio galaxies. Unlike previous studies, we compare radio 
I ' luminosity functions found with both the observed and simulated data to explore the 

, best-fitting fundamental source parameters. The new Monte Carlo method we present 

• here allows us to: (i) set better limits on the predicted fundamental parameters of 
' which confidence intervals estimated over broad ranges are presented, and (ii) gener- 

0^ , ate the most plausible underlying parent populations of these radio sources. Moreover, 

■ as has not been done before, we allow the source physical properties (kinetic lumi- 

| nosities, lifetimes and central densities) to co-evolve with redshift, and we find that 

. all the investigated parameters most likely undergo cosmological evolution. Strikingly, 

CNl ■ we find that the break in the kinetic luminosity function must undergo redshift evo- 

| lution of at least (1 + z) 3 . The fundamental parameters are strongly degenerate, and 

independent constraints are necessary to draw more precise conclusions. We use the 
• i-H ' estimated kinetic luminosity functions to set constraints on the duty cycles of these 

, powerful radio sources. A comparison of the duty cycles of powerful FR lis with those 

determined from radiative luminosities of AGN of comparable black hole mass sug- 
gests a transition in behaviour from high to low redshifts, corresponding to either a 
drop in the typical black hole mass of powerful FR lis at low redshifts, or a transition 
to a kinctically-dominated, radiatively-incfficient FR II population. 

Key words: galaxies: active - galaxies: radio galaxies, quasars - galaxies: evolution 
- galaxies: jets - galaxies: luminosity functions - methods: numerical, statistical. 



1 INTRODUCTION 

Radio galaxies and radio loud quasars are believed to 
have a significant impact on the evolving Univer s e and 
its large scale structu r e (iGopal-Krishna fc Wiital l200ll : 
iRawlings fc Jarvisi 120041 : ISilkl 120051). They are often found 
in galaxy groups and clusters jLongair &; Seldnerl fl979l; 
Hill fc Lillvil99ll:lAllington-Smith et al.lll993l ; lDeltorn et all 
19971 : IZirbell 119971 : iBelsole et a.1.1 |2007| V Since their jets in- 
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ject a significant amount of energy into the surround- 
ing medium, stored in the so-called radio lobes, they 
can provide useful information in the study of the den- 
sity and evolution of the intergalactic and intracluster 
medium. The jet activity is also believed to regulate 



the growth of massive galaxies dBohringer et al.l 



McNamara et a l. 2000; Scannapieco. Silk fc Bouwensl 



1993; 



2005; 



Croton et alj l200tj ; IShabala. Kavirai fc Silkl 120111 ). There- 
fore the estimation of the kinetic power and specification of 
the environments of these powerful radio sources at various 
cosmological epochs is an important task. The intracluster 
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medium can be e xplored with the use of X-ray observation s 
up to z ~ 1 (e.g. iRosati et all 12004 iPoggianti et al.ll2010l) . 
On the ot her hand radio galaxies are found at z ~ 3 and be- 



yond (e.g. Carilli 1 Owen fc Ha rris 1994i; lvan Oiik et al.ll 19961 ; 



Ijarvis et ah 20091 ) and hence they can provide a valuable 



insight into the high-redshift Universe. Moreover, complete 
catalogues of radio galaxies and quasars indicate that these 
sources were of higher number densit y during the so-called 
quasar era i|Dunlop fc Peacock! 1 990), and since they can be 
observed at high redshifts the knowledge of the fundamental 
aspects of these sources as a population is of importance for 
cosmological studies. 

iFanaroff fc Rilevl { 1974) divided the extragalactic 
large scale radio sources into two main classes, namely 
low luminosity FR I and more luminous FR II types. 
Sources of the two classes are morphologically different, 
FR Is are core-jet bright, edge-darkened objects and 
contain presumably turbulent jets, while FR lis are limb- 
brightened, often symmetrical objects with well defined 
features such as jets, hotspots and radio lobes; the latter 
are often referred to as classical double radio sources. 
The question of whether they originate from the same 
parent population but are dependent on the environment 
they reside in, or wheth er they are intrinsically differ- 
ent, is still open (e.g. iBaum. Heckman fc van Breugell 



19921; iBaum. Zirbel fc O'Deal 
Gopal-Krishna fc Wiital |2000|: 



1995: 



Meier et 



Kaiser & Best 



ahf 



Kawakatu. Kino fc Naga3l2009l ; IWang et al.ll201ll ) 



1997 



2007 



Because of their relatively simple structure (contrary 
to the turbulent nature of FR Is), some sophisticated 
semi-analytical models of FR lis' growth have been de- 
veloped. These models can predict source observables, 
that is radio lobe luminosity (L v , where v is the ob- 
served frequency) and linear source size (D), from under- 
lying physical properties such as kinetic luminosity (Q), 
source age (t), and the density of the intergalactic medium 
(p). This has led to a number of studies that investi- 
gated whole po pulations o f these powerful sources and 
their evolution llDalvl [l9 95; Blundcll, Rawlings 'fc Willottl 
19991; iKaiser fc Alexander]|l999l ; iBarai fc Wiital l2006t 120071 ; 
Wang fc Ka iser 2008). The basic methodology of these stud- 



ies was to construct virtual populations of FR II sources 
choosing underlying properties for the population and run- 
ning them through semi-analytical models to predict the 
distribution of source observables, the linear size and ra- 
dio lobe luminosity. The virtual populations were further 
compared to observed data from complete samples of ra- 
dio sources. However, those studies focused on the com- 
monly used radio power - linear size ( P v — D) distribu- 
tion diagram introduced by [Shklovskii (|l963T ). One of the 
major problems with such an analysis is that P u — D di- 
agram is difficult to interpret in terms of most commonly 
used techniques to study source populations such as the 
luminosity function. Moreover, many studies (e. g. Kapahil 



1 19891 ; ISingall 1 1993! ; iBlundell. Rawlings fc WillottlTl999D fo- 
cused predominantly on the relationship between the ob- 
servables and their trends with cosmological epochs. These 
observables are determined by the fundamental source prop- 
erties, and hence they carry convolved effects of possible 
cosmological evolution of the underlying physical proper- 
ties, as well as the influence of observational biases. Al- 
though, there have been attempts to investigate the funda- 




Figure 1. Radio luminosity vs redshift plane for the 3CRR 
(crosses) and BRL (diamonds) samples used in this work (sources 
of FR II morphology with fluxes above the limit S = 10.9 Jy at 
178 MHz (solid line) and with z < 2.0 are shown). 



mental source parameter space (IKaiser fc Alexander! [l999l : 



iMachalski. Chvzv fc Jamrozy||2004al ; IWang fc KaiseiTOOct l I, 
sometimes very strong assumptions have been adopted. 

In this paper we use the standard Monte Carlo method, 
as summarized above, to generate virtual populations of 
FR II sources. However, contrary to previous studies, we 
use distributions of radio lobe luminosity that can easily be 
transformed into radio luminosity functions, rather than the 
P„ — D diagram consistently used in previous semi-analytical 
population studies, to compare the simulated population 
and the observed data. Also, we attempt to make as few as- 
sumptions as possible about the underlying physical source 
properties to obtain more general results (although we are 
still restricted by computing time). To do so we repeat the 
Monte Carlo simulation multiple times following grid min- 
imisation that searches broad ranges of the possible source 
underlying properties. This allows us to create confidence in- 
tervals of the estimated parameters. In addition, as has not 
been done before, we allow for co-evolution of the physical 
source properties; this will enable us to determine the dom- 
inant type of evolution, if any. Further, we investigate the 
influence of various assumptions, such as the density profile 
characteristics or the jet particle content, on the best-fitting 
properties of source populations. The method we develop 
here allows us to constrain the parent population of the 
sources while minimising the inevitable effects of the radio 
surveys' flux limits. 

The paper is structured as follows: the observed samples 
are presented and discussed in !J21 and in Sj3]we summarize 
the theoretical models used. The construction of the simu- 
lated samples and statistical methods used are discussed in 
detail in Sj^l The results are presented and discussed in 33 
and ^respectively, and a summary is given in ij7] 



2 OBSERVATIONAL DATA SETS 

We use complete flux-limited radio samples which contain 
all extragalactic radio sources in a given sky area and above 
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Table 1. Demography of the observational data from the 3CRR 
and BRL radio samples. Quoted numbers of sources valid after 
certain criteria are met (SnsMHz > 10.9 Jy, D > 10 kpc). FR I 
type sources are listed for reference only (for details see ij2]l. 



28 30 
Log 10 (Li78MHz/WHz) 



32 



Figure 2. The observed radio luminosity functions of the anal- 
ysed data (3CRR and BRL catalogues) for each of the considered 
redshift ranges, where Z\ is drawn in black, 22 in blue, and 23 
in green. Each of the redshift bins is further divided in size bins, 
where solid lines (filled circles) correspond to smallest sources, 
dashed (open triangles) to medium size sources, and dotted (open 
square lozenges) to giant sources (see J2]for the exact values). 



the sensitivity limit specified for each survey. Currently, we 
concentrate only on sources of the FR II morphology due 
to the availability of semi-analytical models of their time 
evolution; however, the approach presented here may be 
extended to cover the whole radio source population once 
theoretical models for FR Is are developed. We use both 
radio galaxies and radio-loud quasars as it is assumed, fol- 
lowing the unification models, that the only difference be- 
tween these t ypes of sources is just their viewing angle (e.g. 
Barthcl 1989). Additionally, due to the sparse population of 
sources at high redshifts, and hence poor representation, we 
limit our analysis to sources with redshifts up to 2 = 2. Due 
to the nature of the theoretical models that assume that 
most of the radio luminosity comes from the radio lobes 
of the sources, one should use samples observed at low ra- 
dio frequencies (~ few x 100 MHz) to avoid compact radio 
emission that dominates at GHz frequencies. There are few 
commonly used low-frequency radio catalogues. In this pa- 
per, we present the work done with two such radio samples, 
namely the well known Third Cambridg e Revised Revised 
Catalogue l|Laing. Riley &: Longair 19831. 3CRR), and the 
complete radio sample of lBest. Rottgering &: Lehnertl (| 19991 . 
hereafter BRL). 

The 3CRR radio sample is very shallow (its flux limit 
equals 10.9 Jy defined at 178 MHz), but covers a large area 
of the sky (4.23 sr); it contains some of the most power- 
ful radio galaxies. 146 sources of FR II morphology are in 
the sample; however, we excluded 3C 231 (M82 ), as well 
as, following iBlundell. Rawlings fc Willottl || 19991 . hereafter 
BRW), two additional sources, 3C 345 and 3C 454.5 (due 
to the Doppler boosting). Of the remaining sources 23 are 
of FR I type, and all are found at very low redshifts, not 
exceeding 2 m 0.25. 

The BRL sample is defined at an observing frequency 
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of 408 MHz with flux limit of 5 Jy. The sample has been 
created to complement other radio samples, and specifically 
with the 3CRR catalogue, as it provides coverage of the en- 
tire sky above —30° in declination. Its flux limit translates to 
~ 9.7 Jy at 178 MHz assuming a typical radio spectral index 
a — O. qj. The sample is 100 per cent spe c troscopically com - 
plete (|Best. Rottgering fc Lehnertl l200d ; [Best et alj |20o3 v 
In total, the sample contains 178 sources of which 124 are 
of FR II morphology. 

The two samples are analysed together and both are 
brought to the common observed frequency of 178 MHz. In 
the case of both samples the measured radio spectral in- 
dex of each source is useojf]. The redshift distribution of all 
sources included in our analysis is presented in Figure [1] 
Further, the population is divided into size and luminosity 
bins in given redshift ranges. Due to the number of avail- 
able sources from the catalogues we are only able to con- 
sider three redshift ranges, 21 ^ 0.3, 0.3 < 22 0.8 and 
0.8 < 23 ^ 2.0, the borders of which are chosen to en- 
sure similar source numbers per redshift bin. Incidentally, 
these redshift ranges span over similar relative light travel 
time intervals, i.e. t& zi = 3.370 Gyr, t& Z2 — 3.361 Gyr and 
tAz 3 = 3.364 Gyr. These subgroups of different redshifts are 
considered separately unless otherwise stated. Sources are 
grouped in size bins in the following manner: in 21 the size 
ranges are 10 — 245 kpc, 245 kpc —1 Mpc and > 1 Mpc, in 22 
these are 10 - 208 kpc, 208 kpc -1 Mpc and > 1 Mpc, and 
in 23 10 - 100 kpc, 100 kpc -1 Mpc and > 1 Mpc. The size 
ranges were chosen to contain similar number of sources in 
the case of the smallest and medium sized sources, while the 
last bin tracks giant radio galaxies more accurately (giant 
radio galaxies are defined as sources of total linear sizes of 
more than 1 Mpc). For each of these size bins a radio lobe 
luminosity distribution is constructed. The distribution is 
constructed at at an observing frequency of 178 MHz from 
logio(£i78A/Hz) = 22.0 to logio(Li78Am z ) = 35.0 in steps of 
Alogio(Z/i78Mflz) = 0.5. Such a division of the population 
allows us to track the relative number of the sources with a 
given size, and later to construct radio luminosity functions 
of the sources at the given redshift ranges and of the given 
sizes. The counts of sources used are summarized in Table [TJ 



The relation between the flux of a source (S v ) at a frequency 
v, and its radio spectral index a, is defined as S v oc u~ a . 
2 The spectral indices of the sources contained in the 3CRR sam- 
ple is measured between 178 and 750 MHz, while in the BRL 
sample between 408 and 1400 MHz. 



4 A.D. Kapiriska, P. Uttley & C.R. Kaiser 



and their observed radio luminosity functions are presented 
in Figure [2] 

As pointed out by BRW the use of only one radio sam- 
ple in any population study may result in inevitable bi- 
ases while interpreting results. This is due to the strong 
luminosity - redshift (L„ — z) dependence (Figure [T]). It 
is important, therefore, to either include a few radio sam- 
ples that will supplement each other in the radio lobe lu- 
minosity - redshift plane, or to be very cautious with in- 
terpreting results. For the extensive discussion on this and 
also frequency related issues we refer our reader to the 
origin al BRW pa per. Commonly used catalogues include 
6CE jEaledll985h and 7CRS (base d on the Seventh Cam- 
bridge Redshift Survey Catalogue , McGilchrist et al. 199d ; 



iGrimes. Rawlings fc Will ott 2004; Wang & Kaiser 200|) be 



sides 3CRR (e.g. BRW, IWillott et al.ll200ll ). The BRL 
sample has also been used along with the other three by 
Wang ~fc Kaiserl l|2008l ). However, here, we do not attempt to 
reconstruct the previous population studies of the [P — D — z] 
plane, but rather analyse binned distributions of the sources' 
radio lobe luminosities, hence we present work on only two 
samples. Because of considerably different flux limits of var- 
ious samples, other ones such as 7CRS will be analysed sep- 
arately to 3CRR and BRL. 

We assume a flat Universe with the Hubble constant 
of Ho = 71 km s" 1 Mpc" 1 , and fi M = 0.7 and fi A = 0.3 
throughout the paper. 



Note, however, that the model depends on the 
combination of (p cto) and not the parameters separately. 
The approximation is valid at distances of at least a few 

core radii. 

Further, as iFalld l|l99ll) and KA show, the source ex- 
pansion problem can be solved purely with a dimensional 
analysis where the source growth depends purely on the jet 
of the source with constant power (Q, hereafter referred to 
as radio source kinematic luminosity), its age (t) and the 
environment in which it self-similarly expands. The linear 
length of the lobe (Di) of the source is hence defined as 



A 



ci 



Q 

8 

poOo 



(2) 



where the ci parameter depends solely on /3 and the source 
aspect ratio Rt (i.e. the ratio of the source length to its 
width, see i]4.2.5|l . The lobe pressure (pi) of a radio source 
of a given age is further found with 



Pi CC { Po al ) V3q2/3 jD (-4-«/3_ 



(3) 



To calculate the radio lobe luminosity the KDA model 
incorporates the energy losses of the relativistic electrons, 
using the adiabatic expansion, synchrotron and inverse 
Compton losses. The radio lobe luminosity at a radio fre- 
quency v (note that L v = 4nP^) of the source becomes 
therefore 



3 THEORETICAL MODELS OF RADIO 
GALAXY AND QUASAR EVOLUTION 

A basic picture of the FR II type radio source growth, 
where the relativistic plasma is ejected in two opposite di- 
rections forming collimated outflows, is nowadays widely 
accepted. These outflows expand and interact with the 
surrounding medium until they terminate in a strong 
shock and create so-called radio lobes where the excess 
transferred energy is st ored. The idea first proposed by 
iBlandford fc Reesj (|l974T ) and IScheuerl (|l974h has been de- 
veloped, over the past 15 years, into a few more sophisti- 
cated models of the evolution of FR II sources. In particular, 
there are three mode ls which have gained the most attention, 
iKaiser fc Alexander! dl997|. part 1, hereafter KA) and subse- 
quentlv lKaiser. Dennett- Thorpe fc Alex ander (1997, part 2, 
hereafter KDA), BRW, and lManolakou fc Kirk! j2002| , here- 
after MK). It is not our intent to compare the existing mod- 
els in this paper; however, we will briefly summarize their 
main features. 

The BRW and MK models follow the dynamical evo- 
lution of radio sources as described by KA, which we will 
focus on. The density profile of the environment into which 
the so urce expands, is approximated by the generalised King 
1 19721 ) profile 



L„ tx Qp 



(m + l)/4 



C (x) Ax, 



(4) 



p = p ( — ) for r > a 



(i) 



where r is the radial distance from the AGN core, and p 
is a constant central density within the core radius a . The 
index /3 is usually constrained by the observations, although 
some restrictions may apply in the case of certain assump- 
tions (e.g. the self-similarity assumption requires /3 < 2, see 



where m denotes the injection spectral index of the en- 
ergy distribution of the relativistic electrons f £|4.2,4p . and 
the term C (x) depends on the energy losses and x = ti/t 
with ti being time at which a particle was injected into the 
radi o lobe. For the full derivat ion and discussion see KDA 
and (|Kaiser fc Besdl2007ll20oi ). 

The BRW and MK models differ from KDA in their 
assumptions about the luminosity evolution of the sources 
determined by the way the relativistic particles are injected 
from the jet to the radio lobe, and the particle transport. In 
particular, BRW argue for the importance of the hotspot re- 
gion pointing out the over-simplicity of the KA/KDA model 
and its assumption of source self-similar expansion. Addi- 
tionally, BRW assume that spectral index of the radio lobes 
at a specific time is driven by breaks in frequency in the 
head region. Further, MK expands the KDA and BRW mod- 
els incorporating the particle transport mechanisms. It is 
worth noting that assumptions of BRW and MK lead to 
much steeper evolutionary tracks of the sources in the radio 
power - linear size (P v — D) plane as compared with the 
predictions of the KDA model. 

More recently iBarai fc Wiital <|2006l . 12007m attempted to 
develop a more accurate model by testing and modifying the 
three leading theoretical models. They report that none of 
the existing models can give fully acceptable fits to all of the 
properties deduced from the radio surveys (especially the 
synchrotron spectral index a); however, according to their 
findings KA/KDA give the best overall results as compared 
to the other two models. In this work we employ the original 
KA/KDA model. 
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4 MONTE CARLO SIMULATIONS OF THE 
COMPLETE VIRTUAL RADIO SAMPLES 

In the traditional and most common way, the radio source 
evolution has been investi gated through the so-called P v — D 
diagram (Shklovskii 1963), which in recent studies has been 
extended into radio lobe spectral index (a) and redshift (z) 
parameter space. Here, we use an alternative approach and 
instead of investigating [P v — D — z — a] space we use radio 
luminosity functions (RLFs) in the analysis. An advantage 
of doing so is the direct link to the results of other source 
population studies which are most commonly expressed in 
terms of luminosity or mass functions. 

Due to degeneracy between the fundamental properties 
it is impossible to infer the source kinetic luminosity, the 
central ambient density and the source age directly from the 
observables, and hence Monte Carlo simulations are neces- 
sary in order to generate simulated RLFs. In this section we 
present our Monte Carlo method ( 34. The parameter dis- 
tribution functions and other fixed assumptions of the model 
for constructing the radio source population are discussed in 
34.21 The statistical methods used to determine confidence 
intervals and goodness-of-fit test are discussed in 34.31 



4.1 Monte Carlo simulation 

To construct a single virtual sample we generate, in total, 
~ few x 10 6 virtual radio sources. The subsequent prescrip- 
tion is followed. 

1. For each source of the virtual sample we assign a set of 
physical parameters summarized in Table [2] Each of the 
parameters is either drawn from its respective distribu- 
tion or is the same for each source; the details of these 
physical 'input' parameters are presented in 34.21 

2. Based on the theoretical model of KA/KDA the source 
is evolved and its linear size D is calculated. 

3. The linear size of the source is being checked for its 
reliability, that is whether the source expansion exceeds 
an assumed fraction of the light speed (maximum allowed 
head advance speed, see 34.2.7p . If it does, the source is 
rejected. Otherwise, the source is accepted and its linear 
size is corrected for the projection angle f 34.2,9[) . 

4. The radio lobe luminosity L v of the source at a fre- 
quency v = 178 MHz and its randomly generated redshift 
is calculated. 

5. The radio luminosity is subsequently randomized by 
adding a Gaussian variable of standard deviation AL orr = 
0.08 x L v to account for any instrumental/systematic er- 
rors. Moreover, since the source is evolved in its rest- 
frame, a correction of (1 + z) a transforming the calcu- 
lated radio luminosity to the common frequency must be 
included. The common frequency of 178 MHz is used. 
These corrections are discussed in detail in 34,2.41 

6. Steps 1-5 are repeated ~ few x 10 6 times. 

7. The radio lobe luminosity histograms for each redshift 
range and linear size bin are generated as discussed in 
Ej2j The binning of the virtual sample histograms matches 
exactly the one used for the real observed data sample. 
However, to ensure that the probability density functions 
of kinetic luminosities (Eqn. [8] and Eqn. [9)l are repre- 
sented by reasonable number of sources at the functions' 



exponentially falling ends, we initially generate the vir- 
tual population assuming that their kinetic luminosities 
are drawn from a uniform distribution, and later a weight- 
ing factor is applied to each source (see 34. 2. 111 . 

8. At this stage the histograms represent the total number 
of sources in the simulated sample since we have not de- 
termined yet the fluxes of the generated sources. In this 
sense, those are the 'true' numbers of the entire source 
population. However, to be able to compare the simulated 
sample to the observed data the selection effect arising 
from limited flux sensitivity must be taken into account. 
We achieve it by using radio luminosity functions. RLFs 
are an ideal tool when dealing with coupled selection ef- 
fects such as limited sensitivity of radio surveys. They are 
defined as number density of sources per unit comoving 
volume per unit luminosity. The initial histogram of the 
entire simulated population is therefore transformed into 
RLF by using the survey volume Kurvey, 

, _ A n Li 

<psim,i — A v > \°) 

4t7T Vsurvcy 

where A is the area of the sky that the survey covers, 
and riL i represents the data counts in luminosity bins Li. 
Further, the RLF is transformed back into a histogram 
of number of expected sources for each linear size and 
luminosity bin using Vn m - a maximum co-moving volume 
in which a source in a given luminosity bin and with the 
flux limit of th e survey would be included in the sample 
!|Schmiddlf96Sh , 

47T 

nL i: S lim = —r(j>Bim,iVlim- (6) 
The two transformations simplify to 

nLi,s lita =n Li — . (7) 

V survey 

Note that flux densities of the simulated radio galaxies 
(S'sim.O are not examined directly, and hence no radio 
sources are rejected based on the limited sensitivity bias 
(i.e. requirement that SWm.i > Sum)- Instead, a correc- 
tion factor (Viim/V^urvcy), which indicates the probability 
of radio source occurrence at a given flux limit, is em- 
ployed. The effect of the limited flux sensitivity is an av- 
erage for each considered bin. To ensure greater accuracy 
we perform the flux correction on much finer bin widths 
than those of the final histograms. In particular, in this 
step we use 20 times finer bins than used in the final his- 
tograms. These are later summed to match the initially 
constructed distributions. 

9. So far, we have not discussed the number of progenitors 
becoming active and turning into FR II sources. To en- 
sure a reasonably good statistics, our simulated sample 
is generated with a much larger number of virtual radio 
sources than the observed data sets contain. Instead of 
introducing a corresponding distribution function we use 
the full simulated sample. However, since both the simu- 
lated and observed samples must be of a similar size, that 
is of a similar number of sources considered, we need to 
renormalize the data counts, which are found with Eqn. [7] 
in the virtual sample. We use the binned maximum like- 
lihood method (MLM. see 34.3.11) tog ether with Brent's 
method |Brentlll973l ; |Press et al.lll992h to do so. The nor- 
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malization is set as a free parameter, but represents an 
average for the considered redshift range. 
10. Finally, the goodness-of-fit test is performed. The statis- 
tics used, normalization of the sample and the goodness- 
of-fit test are discussed in i)4.3l 



4.2 Monte Carlo simulation input parameters 

4-2.1 Kinetic luminosity distribution 

The kinetic luminosity of each generated source is drawn 
randomly from a distribution function that acts as a proba- 
bility density function. The form of the distribution function 
of sources' kinetic luminosities is not known, and various 
forms have been assumed in previous works ranging from a 
simple uniform distribution between mini mum (Q m j n ) and 
maxi mum (Q max ) kinetic luminosity (e.g. IWang fc Kaiserl 
2008), through power-law scaling as Q x (e.g. BRW) to more 
complex functions as used by Willott et al. (2Q0ll, discussed 
below). 

We consider two models for the initial distribution func- 
tions of the kinetic luminosities. In Model S the kinetic lu- 
minosity distribution function is assumed to be m odelled by 
the so-called Schechter function (|Schechtedll976h of a form 



exp 



(8) 



where the slope of the function for the kinetic luminosity 
values below the kinetic luminosity break (Qb) is described 
by the exponent a s , and drops exponentially for higher Q. %p* 
is a normalization constant, which in our case is neglected at 
this stage as Eqn.[8]is used as a prob a bility density function. 

Model W follows Iwillott et al.l (|200ll ). who introduce 
modified version of the above Schechter function of a form 



ip(Q) dQ = V* 



Q_ 

Qb 



exp 



dQ, 



(9) 



and where the exponent a s describes the function for kinetic 
luminosities higher than Qb, while for lower kinetic lumi- 
nosities the function drops exponential ly, and ifr* is used as 
in the Model S. In their original paper IWillott et all (|200ll ) 
use a combination of Eqn. [8] and Eqn. [9] to describe the 
whole population that consists of high and low radio lumi- 
nosity sources. Eqn. [9] was introduced to specifically model 
the high radio luminosity subpopulation of radio sources. 
Since we do not consider FR I type so urces in our study w e 
will not follow the original method of I Willott et all l|200ll) . 
and Eqn. [9] alone is used to describe the considered subpop- 
ulation of powerful radio sources. 

Both functions describe the distribution of kinetic lumi- 
nosities between Q m i n and Q ma x, which are set in such a way 
that that contribution to the RLFs of sources with kinetic 
luminosities outside this range is negligible. To ensure these 
probability density functions are represented by reasonable 
number of sources, what will assure good statistics at the 
functions exponential ends, the kinetic luminosities are ini- 
tially generated from a uniform distribution between Q m i n 
and Qmax- Each of these kinetic luminosities is assigned a 
probability of its occurrence (a weighting factor) according 
to the considered probability function (Eqn. [8] and Eqn. [9]). 
This weighting factor is applied to each source while con- 
structing histograms at a later stage. 



4-2.2 Ambient density distribution 

The central density value (p D ) for each generated source is 
randomly drawn from a log-normal distribution with the 
mean value log 10 (p m ) and standard deviation of ci ogl0 ( Po ) = 
0.15. The standard deviation is not introduced as a free pa- 
rameter; however, we test how strong an effect it has on the 
results (see discussion in i)5.4|l . 

Since the type of the surrounding environment (clus- 
ters of galaxies or field galaxies) is tightly linked to the core 
radius of the source, these should never be discussed sep- 
arately. Here, however, we set one value of a for all the 
simulated sources (2 kpc). The choice of a may determine 
the most likely environments found in our simulation; this 
issue is discussed in detail in ij5]and €16.21 The exponent /3 of 
the density profile (Eqn. [1} is randomly chosen from a uni- 
form distribution between j3 m - lYl = 1.2 and /3 m ax = 2.0. Al- 
though ma ny authors u se a co nstant value of P (e. g. /3 = 1.5 
is used bv lDaly|[l995l . BRW. IWillott et al.lll999l ; /3 = 1.9 
and p = 2.0 by KDA, I Wang fc Kaiserl |200S|) ."some obser- 
vational evidence suggest that the parameter may vary be- 
tween sources (e.g. Alshino et al. 2010). 



4-2.3 Age distribution of the simulated sources 

It is assumed that radio sources live up to a certain max- 
imum age (i m ax) after which they instantly die. Although 
this may seem to be an oversimplification as any relic radio 
galaxies are going to be neglected, from previous work of 
KDA and BRW for instance, one can see that the luminos- 
ity tracks of radio sources steepen rapidly at the late stages 
of the source's life, quickly dropping below detectable flux 
levels. 

Here the age of a source is randomly drawn from a 
uniform distribution between t — and imax. The expec- 
tation that all radio sources have the same maximum age 
seems to be unrealistic, hence we introduce a log-normal 
spread of the maximum ages around the mean i ma x m , i.e. 
fiogioOmax) = 0.05. We investigate the effect of the width of 
this distribution in £15.4.21 



4-2.4 Injection and radio spectral indices 

The energy distribution of the relativistic electrons ini- 
tially follows a power-law relation with exponent m, 
N(E)dE oc E~ m dE. The exponent m assigned to a source 
is drawn from a uniform distribution with the minimum and 
maximum value allowed by the theoretical model of source 
growth used in this study, that is m m i n = 2 and m max = 3. 
However, theoretical studies of the particle acceleration in 
the relat ivistic shocks suggest a universal va lue (m ~ 2.2 — 



2.4, e.g. iKirk et al l [2000l ; Ispitkovskvl 120081 ). while studies 



based on observations of gamma-ray burst aft erglows favour 
a Ga ussian distribution of t his p arameter (|Curran et al.l 
120101) . Also, iMachalski et alj (|2007h conclude that the dis- 
tribution of m is rather na rrow with m € [2.0, 2 . 4] acc ord- 
ing to their modelling, and iMeli fc Mastichiadisl (|2008l ) re- 
port on good fits of the spectra by a power-law with index 
m =■ 2.0 - 2.3. 

Moreover, the power-law of the relativistic electron en- 
ergy distribution extends between 7 m i n and 7max, i.e. be- 
tween the Lorentz factors of the least and most energetic 
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Table 2. Overview of a source physical parameters and default distributions from which they are drawn. Details on the assumed 
distribution or value of the respective parameters are discussed in sections as given in column 4. 



Parameter 


Assumed 
distribution 


Description 


Discussed in 






Physical Parameters 
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34 2 21 


il 




V 1 J 1 V 1 f I ' 1 1 Ho 


34 2 21 


8 


distribution 


power law index of the radial density distribution 


34 2 21 


t 


distribution 


source current age 


34.2.3| 


t-inax 


distribution 


source maximum age 




a 


calculated value 


radio spectral index 


34.2.4| 


m 


distribution 


power-law exponent of the relativistic particles' energy distribution 


34.2.4| 


Tmin 


value 


minimum Lorentz factor of relativistic particles 


3|4.2.4| 


7max 


value 


maximum Lorentz factor of relativistic particles 


34.2.4| 


Rt 


distribution 


aspect ratio 


34.2.5| 


■8 


distribution 


projection angle 


34.2.9| 


r x 


value 


adiabatic index of the IGM 


34.2.8| 


r c 


value 


adiabatic index of the radio lobes 


34.2.8| 


r b 


value 


adiabatic index of the magnetic field energy density 


34.2.8| 


k> 


value 


ratio of thermal to electron energy densities in the jet 


34.2.6| 


V max 


value 


maximum allowed head advance speed 


34.2.71 



Distributions' Parameters 



Qb 
"q 

pm 

n T 

^max n 

nt 



value 
value 
value 
value 
value 
value 
value 
value 
value 



kinetic luminosity break 34.2.11 

exponent of the kinetic luminosity distribution 34.2.11 

strength of the kinetic luminosity break redshift evolution 35.2.21 

mean of log-normal distribution of radio sources' central densities 34.2.21 

standard deviation of log-normal distribution of radio sources' central densities 34.2.21 

strength of the central density redshift evolution 35.2.21 

mean of log-normal distribution of radio sources' maximum ages 34.2.31 

standard deviation of log-normal distribution of radio sources' maximum ages 34.2.31 

strength of the maximum source's lifetime redshift evolution 35.2.21 



electrons. We assume 7 m i n = 1 and 7 ma x = 10 10 . KDA 
stresses that 7 m i n <C 7max- We used much higher max- 
imum Lorentz factor than KDA who used 7max = 10 5 . 
Lo rentz factors of 7 ma x ~ 10 5 — 10 6 have also been m e asure d 
bv lMeisenheimer et alj l|l989t ) and lMassaro fc Aielld (|201ll) . 
while BRW favour 7 ma x — 10 14 . We decided to use the in- 
terme diate value as a default one. Further, iBlundell et al.l 
(2006) report a new estimate for the low-energy cut-off of 
the energy distribution of relativistic electrons in FR II type 
sources. Based on the observations of Cygnus A they esti- 
mate a rather high value of 7 m i n = 10 4 . On the other hand, 
invest igations of hot spots done by Meisenh eimer et al.l 
1 19891 ) suggest that the minimum Lorentz factors, below 
which the synchrotron losses are unimportant, are typically 
7min ~ 10 2 , but may reach values of 1 — 10 4 . Values of 10 2 for 
the minimum energy cut-o ff are also supported by, fo r exam- 
ple, iBarai fc Wiital (|2006n and iMeli fc Mastichiadisl (|2008h . 
Since there is a large discrepancy, especially in the estima- 
tions of the Lorentz factors, we investigate the possible ef- 
fects of these different assumptions on the simulated source 
populations in i|5.4.4l 

To find the radio spectral index a one needs at least 
two data points, that is two radio lobe luminosities of the 
same source measured at two different frequencies. A sim- 



ple power-law is employed between L UlTS and Z/„ x to find 
a, which is finally used to correct the L V17S ; the radio lobe 
luminosity estimated with the model, L V17S , is in the source 
rest-frame and one needs to convert it to the common fre- 
quency using (1 + z) a . For consistency and to mimic be- 
ha viour of the 3CRR sample we chose v x = 750 MHz (af- 
ter |Kajse^^^Mexand«i][l999). We n ote here, however, that 
BRW. lJarvis fc Rawlingj (|2000l ) and lJarvis etall (|200ll ) ar- 
gue that this most commonly used relation (i.e. simple 
power-law) may be too simplistic, and instead, curved ra- 
dio spectra should be considered. 

4-2.5 Aspect ratios (Rt) and self-similarity 

One of the consequences of the KA dynamical model of the 
source expansion is its self-similar growth. Instead of mea- 
suring or introducing the jet opening angle, the aspect ratio 
(Rt) is used. In the KA model the aspect ratio is linked 
to the jet opening angle (6) by -Rt oc 1/9 based on the as- 
sumption that the lobe pressure is balanced by the external 
gas ram pressure. _Rt stays constant through the source life- 
time. The BRW model, on the other hand, assumes that Rt 
increases with the source growth, and their results suggest 
some dependence on the source kinetic luminosity and/or 
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its age, while iMachalski. Chyzv fc Jamrozvl (|2004al lbh sug- 
gest th^t_th^_sdf;simija^ not h old for old sources. 

iMullin. Rilev fc Hardcastld (|2008l ) present a detailed in- 
vestigation of aspect ratios of 3C sources of z < 1.0. The ob- 
served Rt values are found to fall in a range of 1 < Rt < 8, 
with a median within (Rt) £ 1.6 — 2.6 (which depends on 
spectral class). Mullin et al. (2008) report that higher as- 
pect ratios seem to occur for larger sources (> 100 kpc) 
which may suggest a non self-similar source growth, or a self- 
similar growth occurring only in the early stages of a source 
life. However, the aspect ratio has a predominant influence 
on the radio source linear size (hence its age). In particular, 
the smaller the aspect ratio gets, the wider the opening an- 
gle of the jet becomes, and thus a higher pressure of the lobe 
is needed to account for the faster sideways expansion of the 
lobe. This in turn will lead to smaller head advance speeds 
of the jet, as well as larger synchrotron losses, and will result 
in smaller linear sizes as compared to the sources with larger 
Rt- This may be one of the reasons why larger sources seem 
to have larger aspect ratios. Similarly, BRW suggested that 
the aspect ratio might be higher for more powerful sources. 
One must notice, however, that for a given age and envi- 
ronment of a radio source higher kinetic energy will lead 
to larger linear size, and, according to our argument above, 
larger aspect ratios may be expected. 

Since the distribution of _Rt is not yet well constrained; 
we decided to assume a uniform distribution of the aspect 
ratio with 6.0 > -Rt > 1.3, which translates to the range of 
jet opening angles of 13.7° < 8 < 63° (values based on 
iLeahv fc Williamslll984l ; lDalvlll995l ; Machalski et al. 2004a; 
Mullin et al. 2008), and to follow the original KDA model 
which implies self-similar growth of a radio source. A dis- 
tribution of aspect ratios allows for a variation in growing 
rates of radio sources and hence is more realistic than a sin- 
gle value. 

4-2.6 Jet particle content 

The ratio (k 1 ) of the energy densities of thermal particles to 
the energy densities of the electrons at the time they are in- 
jected into the cocoon is assumed to be k! — 0. Note that this 
definition of k! differs slightly from typical assumptions (k) 
in such a way that relativistic and non-relativistic electrons 
are already included even if k' — (this already implies that 
the typically assumed ratio is k ^ 0). 

There is much debate on the particle content of the 
radio galaxy jets. Some argue that the FR II jets are 
lightweight (electron-positron dominated), while FR I jets 
are heavy, that is they may possess a significant pro- 



ton content ( e.g. ICelotti fc Fabianlll993];IWardle et allll998l; 



Kinq fc Takaharall2004l ; lpunn fc Fabiarj2004ICroston et ail 
20051 . 1200881 ). Moreover. Illardcastle fc Crostonl (|2010h based 



on the investigation of Cygnus A report that FR II radio 
galaxies may be more proton dominated (fe ~ 1 — 4) if they 
reside in very rich environments. KDA showed that heavy 
FR II jets will require significantly higher jet powers to re- 
produce the same radio luminosity as the electron-positron 
dominated jets, hence they concluded that k' must be close 
to 0. Initially, in our simulation we followed the conservative 
KDA assumption, but the effect of changing k' (to allow a 
proton content in the jets) on the whole population of sim- 
ulated sources is investigated in £15.4.51 




Figure 3. The goodness-of-fit test. The synthetic log-likelihood 
distribution of the 2 X 10 5 synthetic data sets generated from 
the model average histogram (consult £|4.3l for term explanation). 
The red arrow indicates the position of the actual data set log- 
likelihood. The dashed line (black) marks the 68.3 per cent con- 
tainment (from the minimum value of — ln(£), note that the distri- 
bution is one-sided), and the dot-dashed one (black) the 95.4 per 
cent containment. The histogram of the synthetic log-likelihood 
distribution of the lowest redshift data of Model S is shown. The 
data is consistent with the model at the 90 per cent confidence 
level (p-value=0.234). 



4-2.7 Maximum head advance 

It is at first assumed that the head advance speed of a source 
("Umax) at the time of observation may not exceed 0.4c, all 
sources that surpass this limit are rejected at the stage of 
generating the population. Such an upper limit is consis- 
tent with the speeds inferred from the synchrotron spec- 
tral ageing of high luminosity double rad io sources (e.g. 
iLiu. Poolev fc Rilevl[l992l ; iBest et al.lll995l ). and supported 
by the estimation of d ynamical ages of FR II sources by 
IMachalski et al.l (|2007l ) who report f max ^ 0.3c. However, 
there have been discussions on the possible overestimation 
of the lobe advance speed upper limits; results obtained 
through the analysis of lobe asymmetries and the steepen- 
ing of radio spectra su ggest that the head advance speeds do 
not ex ceed 0.15c (e. g. Ashkarian fc Longald Eoool ) . or even 
0.05c (|Scheuerlll995r i. We discuss the effect of different as- 
sumptions of Umax on the results in £15.4.31 



4-2.8 Adiabatic indices of radio lobes, magnetic field, and 
external medium 

After KA/KDA, we assume the adiabatic indices of the radio 
lobes, magnetic fields and the external medium to be Ti = 
Tb = 4/3 and T x = 5/3 respectively, that is we assume 
a non-relativistic equation of state for the closest external 
medium of the radio source, and relativistic particles inside 
the radio lobes. 



4-2.9 Projection effects 

The linear size calculated based on the randomly chosen and 
set physical parameters of each source in the generated pop- 
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Table 3. Searched ranges and steps of the distribution parameters in the grid minimisation. 
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ulation (see steps 2. and 3. in EJ47TJ is the source true size, 
and is further randomly oriented in the sky as it is observed 
at some projection angle It is assumed that the projection 
angle is distributed uniformly in [1 — cos(#)] plane. There- 
fore, the size becomes -D pro j = Z3trueSin(#). For simplicity 
we will refer to this projected size as D, and the true size is 
not considered from now on. 

Also, note that Eqn.[T]is an approximation and is valid 
for distances r greater than few core radii, hence we do not 
consider sources smaller than 10 kpc in our samples, either 
the simulated or the observed oneaj. 



4-2.10 Redshift distribution 

It is assumed that sources are uniformly distributed in space, 
between minimum (Knin) and maximum (Vmax) comoving 
volume depending on the considered redshift range. From 
this a radio source redshift (z) is found. Redshifts from z — 
to 2 = 2 in steps of Az = 0.001 are considered in this work. 
Note that evolution of source number density is not initially 
modelled in the construction of the simulated population; 
however we apply a normalization (scaling), which is a free 
parameter in each redshift range (step 9. in i)4.1l t)4.3.1[) . 
before the maximum likelihood method and goodness-of-fit 
test are performed. Hence we obtain an average in radio 
source number density for each z bin considered. 



4.3 Confidence intervals and Goodness-of-Fit test 

4-3.1 Binned Maximum Likelihood Method 

We use the binned maxim um likelihoo d method (MLM) in 
the statistical tests (e.g. ICashl Il979l ). A log-likelihood is 



3 9 sources in total, resulting from this requirement, were ex- 
cluded from the observed samples. 



found in fitting each of the simulated samples to the ob- 
served data according to 

hxC = E^r [sciln g ) - | - In(x 4 !)] (10) 

where N is the number of bins, each with an expected num- 
ber of sources yi and with the observed number of counts 
Xi, and n is the normalization constant discussed in step 9. 
in $Q] 

4-3.2 Goodness-of-Fit 

To perform the goodness-of-fit (GoF) test we use the Monte 
Carlo analysis. The final histogram of the simulated radio 
sources is in fact the mean of all possible realizations of the 
underlying 'true' population, which we call the 'model av- 
erage histogram'. The Monte Carlo procedure undertaken 
here is to randomize each bin of the model average his- 
togram within the Poisson regime to obtain a 'synthetic 
data set', a single realization of the 'true population'. Fur- 
ther, we apply the binned MLM method to find the log- 
likelihood between the newly created synthetic data set and 
the model average histogram. We repeat the procedure mul- 
tiple times (~ 2x 10 5 ) to obtain a distribution of the possible 
log-likelihoods of the true population, which we refer to as 
the 'synthetic log-likelihood distribution'. Finally, the log- 
likelihood of the observed data histogram is compared to the 
generated synthetic log-likelihood distribution. The la level 
of consistency requires the actual data set log-likelihood to 
be placed within the 68.3 percentile of the log-likelihood 
distribution, i.e. between and 68.3 per cent since our dis- 
tribution is one-sided (Figure [3} • We quote the GoF test 
results in terms of p— value, which is the probability that 
the test statistic is at least as extreme as the one observed 
and assuming the null hypothesis is true, la is equal to 0.317 
in terms of p— value. The null hypothesis is rejected if the 
p— value is less than the significance level a s i g ; we assume 
a s i g = 0.1 = 10 per cent. The choice of the number of syn- 
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thetic histograms allows us to check the agreement up to 5a 
level (p — value = 6 x 10~ 7 ). It is important, however, to 
remember that the confidence levels obtained in such a way 
are nominal only as the synthetic log-likelihood distribution 
is not Gaussian in our case. 

Although this technique tests whether a model fits data 
well, it is incapable of distinguishing the best-fitting models. 
To compare these best-fitting models one needs to perform 
the so-called likelihood ratio test the statistic of which is 

d=-2[ log(£(# m0 )) - log(£(# m i)) ] (11) 

and where £ denotes log-likelihood, H m o is the null model 
and H m i an alternative one. The p- value of the d statistics 
may be obtained to find which model should be preferred. 
The models must be nested, however, and if this requirement 
is not fulfilled, the likelihood ratio test is invalid. For the 
models -ff m o and fl^i to be nested it is necessary that if m i 
contains the same parameters as H m o, and has at least one 
additional parameter. 

4-3.3 Confidence intervals 

The source physical parameters that are the main focus in 
this study are t maXm ( E|4.2.3|) . p m ( i|4.2.2[) . and Qb and a s 
(i j4.2.1l) . To obtain the confidence intervals of these param- 
eters we perform grid minimisation searchin g thei r broa d 
ranges (Table [3]), and follow the method of ICashl l|l979l) . 
For each parameter set in our searched grid there is a cor- 
responding log-likelihood. From these, the global £ max is 
found, indicating the best fit to the ob s erved data (Eqri ll0[) . 
Furthermore, as pointed out by I Cash! |l979), one may find 
the difference between £ max and £ for all the other sets of 
parameters in the evaluated grid, the so-called AC statistic, 
which is defined as 

AC = [-2 ln(£ max ) + 2 ln(£)]. (12) 

£ is the log-likelihood of the sub-grid, which is extracted 
from the global grid by setting each point of a parameter 
which we are focused on as non- varying and the correspond- 
ing set of log-likelihoods of these points are listed based on 
the all the other parameters that vary. 

The AC statistic is distributed as % 2 with n£ max — nc 
degrees of freedom (dof), where nc max denotes degrees of 
freedom associated with the global £ max and nc denotes 
the dof of restricted sub-grid (dof = 2 in our plots). The 
confidence intervals are, therefore, defined such that con- 
tours encircle parameters for which their log-likelihood is 
above a certain value (£ > £o), and levels of \ 2 distribution 
may be used. Note that this may lead to disconnected re- 
gions in the case the likelihood function is highly irregular. 
For an in-depth description of constructing co ntour plots in 
case s such as this one we refer the rea der to ICashl (|l976T ) 
and lLampton. Margon fc Bowverl l| 19761 ). 



5 RESULTS 

In this section we present and discuss the results of our 
Monte Carlo simulations. The fitted histograms and the ra- 
dio luminosity functions of the analysed observed and simu- 
lated samples are presented in i]5.1l The evidence for the cos- 
mological evolution of the physical parameters is discussed 



in i]5.2l An important note on the parameter degeneracy is 
discussed in i|5.3l Finally, ^5.41 contains discussion on the 
possible effects of the model assumptions on the results. 

5.1 Radio luminosity functions and fitted data 
counts 

The best-fitting model histograms are presented in Figure [4] 
The corresponding RLFs are shown in Figure [5] The quoted 
uncertainties on the RLFs are the Poissonian errors only. 
The simulated RLFs appear to be consistent with the ob- 
served data. Moreover, we were able to reconstruct the lumi- 
nosity distributions of the observed samples for their given 
linear size subsamples as well as the redshift ranges. The 
match in the linear size distribution is of particular inter- 
est, as it has been noted previously that not all previous 
work succeeded in reconstructing the linear size distribu- 
tions an d often too many large so urces have been created 
(see e.g. lBarai fc Wiitalliooel . 120071 ) . It may be possible that 
such an effect might not be due to the specific models used, 
but rather is due to restricting the model parameters, which 
in our case were kept free. We have not tested analytical 
models other than KA/KDA in this work. 

5.2 Evidence for the cosmological evolution of the 
intrinsic and extrinsic source parameters 

In the simplest case we have tested, it was assumed that the 
parameters that are our main focus (that is t maXm ,Pm,QB 
and a s ) do not evolve with redshift within each redshift bin. 
We analyse the redshift ranges independently, hence sepa- 
rate sets of best-fitting parameters for each redshift bin are 
found; if there is no redshift evolution of the intrinsic and ex- 
trinsic parameters, the results for each redshift range should 
not be significantly different. These are referred to as the 
independent-redshift fits and are presented in i]5.2.1l Sub- 
sequently, we have attempted to investigate the strength of 
the cosmological evolution of the source physical parameters 
by allowing for continuous redshift evolution of t maXm , p m 
and Qb within each redshift bin, and fitting all the redshift 
bins simultaneously using the same parameters to describe 
the evolution. We refer to this case as the combined-redshift 
fits, and the results are presented in i]5.2.2l 

5.2.1 Independent-z fits 

The confidence intervals of the best-fitting parameters for 
the simple Model S in the respective redshift bins are shown 
in Figure [6] - Figure [9] and similarly for the Model W in 
Figure [TD]- Figure [T3] The best fits and the GoF test results 
are listed in Table [4] although we once again stress that 
the best-fitting parameters maximize the likelihood but are 
not the standard means of their respective distributions and 
the inspection of the associated confidence intervals is very 
important. 

The confidence intervals span a wide range of the pos- 
sible values for each of the parameters resulting in relatively 
large uncertainty in the best fit. One must bear in mind that 
these confidence intervals are based on the underlying five 
degrees of freedom allowed in the simulation, and any visible 
degeneracy or shift may be influenced by a change of other 
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parameters. Indeed, some degeneracies are seen, and are es- 
pecially strong in the case of the maximum source age and 
the corresponding central density; however, the degeneracies 
seem to be also present for the Qb — tmax m and Qb — p m 
results. The degeneracy issue is discussed separately in £|5 . 31 
The confidence intervals do not converge for the high 
kinetic luminosity break values (> 10 42 W) in the case 
of Model S, and for the low kinetic luminosity breaks (< 
10 37 W) in the case of Model W. It is accompanied by neg- 
ative slopes of approximately 0.5 — 1.5 in both cases. At 
this stage the function becomes a power-law. Despite the 
fact that the kinetic luminosity break indicates values of 
> 10 42 W, the number density of jets with kinetic lumi- 
nosities close to this break is so low that the actual Q max 
contributing to the simulated population is ~ 10 41 W. Sim- 
ilarly, in the Model W the kinetic luminosity break shifts 
to very low values and only the high kinetic luminosity end 
of the distribution, which at this point is approximated by 
a power-law, contributes to the observed data counts. In 
previous studies, authors have assumed such power-law dis- 
tributions. Assuming that 



iP{Q)AQ = Q- x dQ, 



(13) 



slopes of i = 2.6 (BRW), x = 3.3 and x — 3.6 
ijBarai fc Wiital 1200611 . x = 16 llKaiser fc Bestll2007l ). and 



x ~ 0.9 - 1.3 jWang fc Kaiserll2008l ) have been found. Here, 
we observe slopes of x ~ 0.5 — 1.5 (Model S and W) in the 
power-law extremum. We can conclude here, therefore, that 
the hypothesis that kinetic luminosities follow an unbroken 
power-law distribution cannot be ruled out at a confidence 
level of more than 95.4 per cent. 

To ease the investigation of the possible cosmological 
evolution of all the searched parameters Figure [9] and Fig- 
ure \T3\ show superimposed 90 per cent intervals of all red- 
shifts. We find that the kinetic luminosity break shifts to 
higher values for higher redshift ranges in the case of both 
Model S and Model W. Furthermore, Figure [5] clearly shows 
that there is no common solution found for the maximum 
source age and the central density for all three redshift 
ranges. Consequently, the hypothesis that there is no evolu- 
tion of tmax m and p m is ruled out with probability of > 99 
per cent. One may notice, however, that it is possible to find 
Qb to be constant for all redshifts (Model S only), but this 
would imply unrealistically strong evolution of p m . To avoid 
such a strong evolution more than one parameter would have 
to undergo an evolution with redshift. As an additional test, 
we have attempted to fit the three redshift bins simultane- 
ously with one set of parameters that would be valid for all 
redshift sub-samples. No satisfying results have been found 
for the latter test case, the consistency of the hypothesis 
with the data was ruled out at the nominal 5<r level. 

5.2.2 Combined-z fits 

We have attempted to investigate the strength of the cos- 
mological evolution of the intrinsic and extrinsic source pa- 
rameters. The results of the independent- z fit suggest evo- 
lution of more than one parameter, and bearing in mind 
the possible degeneracies (see discussion in i|5.3[) . we tried 
to restrict as few parameters as possible. We allowed simul- 
taneous evolution of the following parameters: (i) the ki- 
netic luminosity break was allowed to evolve with redshift 



as Qb(z) — Qb(z — 0)(1 + z)™ q , (ii) the central density of 
the environment as p m (z) = Pm(z = 0)(1 + z) n ' , and (iii) 
the maximum age of the sources was assumed to undergo 
evolution according to t maXm (z) = t mB0Cm (z = 0)(1 + z) nt . 
However, due to limited computing time and power, we were 
forced to compromise on the resolution of the grid of the 
parameters searched over. For instance, the n q ,n r and nt 
exponents were searched in steps of 0.5 only, and QB,p m 
and t m ax m in steps of 0.3 in a logarithmic scale. 

The results are displayed in Figure [15] (Model S) and 
Figure [16] (Model W). The best agreement between the sim- 
ulated populations and the observed samples is listed in Ta- 
ble [4] As discussed in 35.31 the degeneracy between param- 
eters means that the n q ,n r and n t cannot be easily con- 
strained, and additional constraints on the parameters are 
necessary. An interesting result emerges for the kinetic lu- 
minosity break cosmological evolution, for which the trend 
with redshift seems to be unavoidable (n q > 3) as our results 
show. We discuss these results in detail in 36] 



5.3 Central density — maximum age degeneracy 

The degeneracy between central density and maximum age 
needs some extra attention as it seems to be often over- 
looked. The calculated radio lobe luminosity of a radio 
source depends on its kinetic luminosity, central density and 
the age of the source. However, one must be aware that re- 
stricting one parameter will yield an adjustment in another. 
In particular, the maximum age of the sources and the cen- 
tral density in which sources evolve seem to be the most 
strongly correlated. From Eqn. [3] -[3] one can deduce that 



L v oc 



2 l B -0) (/J fl£)3T5=?Tt sis 



6-4m-ff(3 + m) 



C{x)d(x) 



Assuming m = 2.5 and /3 = 1.5 (most commonly assumed 
values) and ignoring for the moment the energy losses term, 
one will see that 



r ^,1.25 1.15,-1.75 

L v oc Q p t 



(15) 



Clearly, if we keep radio lobe luminosity, kinetic luminosity 
and all the other values the same for a source, and change 
only its central density and the age, the two latter parame- 
ters will compensate for each other. For example, a change 
in the central density by a factor of 10 (reaching less dense 
environments) will yield a change in source age by a factor of 
7 (implying younger ages) to maintain the same radio lobe 
and kinetic luminosities. Of course this is an approximation 
only as such a change will also yield changes in the source 
linear size and in the loss processes term C(x) of Eqn. [3] 
since both depend on the stage of the source life. Note that 
Eqn. 1151 cannot be used as an approximation for simple L v 
calculation because the energy losses cannot be neglected. 
Eqn. [15] is used here solely for descriptive purposes. 

One may also argue that the kinetic luminosity may 
compensate for the change of the source age instead of the 
ambient density (Eqn ll4[) . Our results suggest, however, that 
the p m — £ maXm degeneracy is dominant. The effect of this 
degeneracy is seen in Figure [B] - Figure 1121 and addition- 
ally in Figure Q2] where results in the Qb — p m plane for 
source populations of different assumed maximum ages are 
displayed. 



Table 4. The best-fitting parameters for the all tested cases of Model S and W for each redshift range. Due to occurring degeneracies (§5.3) one should 
always consult the corresponding confidence intervals (sec Figure 6 — Figure 13). The following standard deviations of logio(p m ) and logio(t maXm ) log-normal 
distributions are used: logio(o> ) = 0.15 and logio (<xt max ) = 0.05. 90 per cent uncertainties are quoted. 



to 

b 



Model 



logio(Q B /W) 



l°glo(Pm/kgm 3 ) 



log 10 (i m ax m /yr) 



Model S 

independent-z fits 



Model S 
combined-2 fits 

Model W 
independent- 2 fits 



Model W 
combined-2 fits 



21 
22 



23 



all 



21 



23 



all 



39.15 
39.00 
39.90 



0.30 
-0.30 
0.60 
0.45 
+0.60 
0.45 



38.0+( < 3 °- 3) * 



38.25t°-« 
40.601°;?° 
41.35l°;f 



37.8 



+«o.3)t 
+«o.3)t 



10.5 



+ (<0.5)t 
-0.5 



0.3 
-0.6 
+ 1.3 
-7.7 



0.6 
-0.9 



0.6 



-«0.2)t 
-0.6 



i-so±S:5g 

7.30l 5 5 j°t 

1+7.401 
-3.95f 



3.50 1 



2.8 



+ «0.2)J 
+ «0.2)t 



-23.4l°;« 



-20.4 
-20.0 



+ 1.8 

2.0 
+2.0 

2.8 



„, „+«0.3)t 
-^• u _(< .3)t 



-23.4l°;g 
-22.6l°j 



-23.0 



+2.2 
0.4 



-23.6 



+ «0.3)t 
-«0.3)t 



n „+«0.5)t 



1.0 



+ «0.5)t 
-«0.5)t 



7.23 
7.83 
7.13 



-0.15 
0.60 
0.75 
+0.60 
1.05 



7 „+«0.3)t 
'■°-«0.3)t 



7.23 
7.08 
6.08 



+0.45 

0.45 
+0.30 
0.30 
+0.75 
0.15 



,+ «0.3)} 
5 -(<0.3)t 



_ 4 Q +«0.5)t 
4 - U -«0.5)t 



-5.0 



+ «0.5)t 
«0.5)t 



p— value 



0.234 
0.639 
0.925 

0.380 



0.251 
0.674 
0.885 

0.557 



to 
?r 



-5 



Notes. The resolution of the results is Alog 10 (p m ) = 0.2, Alog 10 (Qe) = 0.15, Alog 10 (t max ) = 0.15, and Aa B = 0.1 for the independent-2 fits, and Aa B = 

0.2, Alog 10 (p m ) = Alog 10 (QB) = log 10 (tmax) = 0.3, and Ant = An r = An q = 0.5 for the combined-2 fits. 

f For errors which may be extending beyond the searched ranges (see Table 3) value up to the range border is quoted. 

X If errors are smaller than their respective resolution, the value of < A is quoted. 
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5.4 Influence of the assumptions on the other 
model parameters 

Although, one can easily estimate how strongly a given phys- 
ical model assumption influences a single source, such pre- 
dictions may not be trivial when considering the whole pop- 
ulation of radio sources. Here, we discuss the possible influ- 
ence of the model assumptions introduced in H4.2I on the 
simulated source populations and our results. The corre- 
sponding confidence intervals for each case are plotted in 
the online Supplementary Material. 



5.4-1 Core radius and the shape of the central density 
distribution 

The central density is tightly linked to the core radius a , 
and, we stress again, the parameter (p af) should not be 
considered separately in terms of p and a . In our simula- 
tion the core radius is kept fixed at a value of 2 kpc, and 
the corresponding values of the mean central density are 
searched through. Any change to the value of a will affect 
the result of p only, where for lower values of the core radius 
the mean central density will compensate shifting towards 
higher values, and for larger a less dense environments will 
be preferred (assuming no change is made to p). 

The initial log-normal distribution of central densities 
is assumed to be quite narrow, with er logl0 ( Po ) = 0.15. We 
tested how much the results will be affected if we allow 
a broader central density spread. As we discussed in £15.31 
the change in density will have very strong effect predomi- 
nantly on the maximum allowed age of the sources. More- 
over, since both radio luminosities and linear sizes depend 
on the density of the medium, allowing for larger standard 
deviation of the p m distribution will affect the distributions 
of kinetic luminosities of the sources and their ages as these 
will have to compensate for the environment change to re- 
produce the observables (as compared to the initial narrow 
cr logl0 ( Po )). Indeed, the results show that the kinetic lumi- 
nosity break shifts to slightly lower values as the central 
density standard deviation gets broader. The mean central 
density and the maximum age of the source of the best fits 
seem to oscillate around similar values for the tested cases, 
where o"i O gi ( Po ) = (delta function), o" logl0 ( Po ) = 0.15, 
°"iogio(po) = °- 50 and °"iogio(Po) = °- 75 I in units of logio(Po)]. 
However, we see a significant change for 0"i O gi O (p o ) = 1-00, 
when the confidence intervals of all parameters become more 
defined, in particular Qb drops to values of 10 38 W and the 
maximum age of the sources oscillates around ~ 1 x 10 s yr 
in the case of Zi, while for zz there are only very specific 
regions of the parameter space allowed. 



5.4-2 Distribution of the maximum source lifetimes 

The initial log-normal distribution of the maximum 
source lifetimes is assumed to be very narrow, with 
""logici (Wx) ~ 0.05. We tested how much the results will 
be affected if we allow for greater variety in the maxi- 
mum lifetimes within one population of the sources. We re- 
run the simulation twice, with o"iogi (t max ) = 0.3, and with 
°"iogio(t m£ , x ) = 0.6. Similarly to the case of the standard devi- 
ation of the central density distribution, we do not observe 



any significant changes in the confidence intervals of any of 
the searched parameters. 

5-4-3 Head advance speeds 

Since there is much discussion on the value of maximum 
head advance speed of the FR II jets (see £|4.2.7[) we re-run 
the simulation with f ma x = 0.05c, u max = 0.15c, and further 
with u m ax = 0.8c, to investigate how strong an effect such 
an assumption may induce on the whole population. We no- 
tice that the lower u max is the more constrained confidence 
intervals become. Although there is no significant difference 
between the confidence intervals, stronger constraints, re- 
sulting from lower allowed w max , indicate that higher cen- 
tral densities, older source ages, and slightly less powerful 
sources (the latter seems to be the case only in the smallest 
redshift range) from the initial broad contours are indeed 
preferred. 

Such results come from the fact that, in the non- 
relativistic case, the energy density in the jet head (hot spot) 
is mhs = Px^advi where p x is the external medium density, 
z^adv is the head advance speed, and f a dv = ^/[l + CPx/pj) 1 ^ 2 ] 
with the jet speed Vj and jet density pj (see Marti et al. 1997 
for a relativistic extension of this calculation which leads to 
very similar results). This implies that u^dv oc p x (see 
IScheuer|[l996h . Thus, the smaller the head advance speed is, 
the higher the ambient densities of the sources must become 
to reproduce their observed linear sizes and radio luminosi- 
ties. If the external density changes then adjustment in the 
tmax and Qb must develop to again ensure that the observed 
L v and D are reproduced. This is indeed what is observed. 

Note, however, that since from the KA model we know 
that the expansion speed decreases with the source age, the 
maximum head advance speed that we refer to here is in 
fact the maximum expansion speed a source may have at 
the time of observation. Since there is a higher probability 
of observing an old source than a young one, these t> max may 
not be the highest possible expansion speeds during a source 
lifetime as at the earlier stages of the source life its early-i/ a dv 
(non observed) may be higher than the decelerated late-f a dv 
(observed, at the later stages of the source life). The latter, 
only, is compared to w ma x. 

5-4-4 Energy distribution of the relativistic electrons: the 
injection index and the Lorentz factors 

Initially, we assumed a uniform distribution of m £ [2, 3] . 
We tested, however, the possible influence of the assumed 
m distribution on our results, and we re-run the simulation 
twice assuming that m is set to a single, universal value 
(m = 2.3), and further that m follows a Gaussian distribu- 
tion with m me an = 2.4 and cr m = 0.3. Note that this nor- 
mal distribution is asymmetric here since it is not allowed 
to extend beyond t he minimum and maximum m values. 
ICurran eT al. (2010) finds m mC a n = 2.4 and a m = 0.6; how- 
ever, since m is restricted by its minimum and maximum 
allowed value in our study, such a distribution would be 
nearly uniform for m m i n = 2 and m max = 3. For this rea- 
son we decided to assume a smaller standard deviation of 
the distribution in our tests. We find no difference between 
results found with either of the Gaussian, the uniform dis- 
tributions or a single value of the injection index. Hence, 
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Figure 4. Histograms of the observed radio sample (black solid line) and simulated ones created with the best-fitting parameters of 
the independent-^ fits of model S (black dashed line) and model W (black dotted line), and the combined- z fits of model S (blue 
dashed line) and model W (orange dotted line), in z\, Z2 and 23. In each redshift range there are three separate size bins as described 
within the subplots. Each redshift bin is simulated independently, while the FR lis linear sizes within each redshift range are simulated 
simultaneously; that is, a good fit to all linear sizes simultaneously at the same redshift is required. The simulated source populations 
created with the best-fitting parameters are consistent with the data at the 90 per cent confidence level based on the AC statistics (for 
exact p— values see Table [4}. 



we conclude that we are not sensitive enough to distinguish 
between the underlying distributions of m at this point. 

Furthermore, since there has been discussion on the low- 
energy cut-off of the energy distribution of relativistic elec- 
trons present in the radio jet we tested the following cases to 
investigate the possible changes to the results: (i) 7 m i n = 1 
and 7 max = 10 5 , (ii) 7 min = 10 2 and 7, nax = 10 10 , and (iii) 
7min = 10 4 and 7 max — 10 10 . All of these were compared to 
our initially assumed case where 7min = 1 and 7 max = 10 10 . 
As expected, there is no difference between the initial as- 
sumptions and case (i), indicating that indeed the exact 
value of the maximum Lorentz factor is not crucial as long 
as 7min >C 7max- In case (iii) we observe a significant change 
of the kinetic luminosity break confidence intervals, where 
for 7 m in = 10 4 the kinetic luminosity break is smaller by 



~ 1.5 decades, as compared to the other cases. Such a dras- 
tic change in confidence intervals in case (ii) is not observed, 
despite 7 m i n being significantly larger than our initial case. 
Here, we conclude that change in the minimum Lorentz fac- 
tor for the whole source population will significantly affect 
their kinetic luminosities only. 

7min, 7max and m directly influence the initial energy 
density distribution of the relativistic particles. The mini- 
mum Lorentz factor indicates how relativistic are the least 
energetic particles. The lower 7 m i n the more cold material is 
included. For cold plasma higher Q is required in order to ob- 
tain the observed radio luminosity because a fraction of the 
kinetic luminosity will be lost for particle acceleration. How- 
ever, if this material is not included, less power is required 
to reproduce the observed L, and hence significantly lower 
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Figure 5. Radio luminosity functions of FR II type sources from 3CRR and BRL catalogues (solid) and the simulated populations 
generated with the best-fitting parameters for each redshift and model tested (independent-z fits of Model S drawn as black dashed line, 
independent-z fits of Model W as black dotted line, combined-z fits of model S as blue dashed line, and combined-z fits of model W as 
orange dotted line). Each redshift bin is simulated independently, while the FR lis linear sizes within each redshift range are simulated 
simultaneously; that is, a good fit to all linear sizes simultaneously at the same redshift is required. The simulated source populations 
created with the best-fitting parameters are consistent with the data at the 90 per cent confidence level based on AC statistics (for exact 
p— values see Table |4j. 



Qb are observed in our results if we allow 7 m i n to be as high 
as 10 4 . The maximum Lorentz factor characterizes the cut 
off at the high energy end of the particle distribution; parti- 
cles can only be accelerated up to a certain frequency above 
which loss rates take over (see also lMarkoff. Falcke fc Fended 
2001). The index m defines the slope of the energy density 
distribution, and hence determines the relative fraction of 
low and high energetic particles, where the flatter the distri- 
bution is the more energetic particles are included. Although 
the change of m may be significant for a single radio source, 
it does not seem to have a dramatic effect on the whole 
population of sources, and it is the low energy cut off that 
influences the average properties of the population. 



5.4-5 Jet particle content 

As discussed in i]4.2.6l there is much debate on the particle 
content of the radio jets. Initially, in our simulation we fol- 
lowed the conservative KDA assumption (k' = 0), but to test 
the influence of this assumption on the whole population of 
sources we assumed the particle content of radio jet to be 
drawn from a uniform distribution in a range k' 6 [0, 10], 
and further also k' £ [0, 100]. As an additional test we also 
investigated the case when k! = 100. A significant change in 
the fitted source parameters were observed only in the two 
latter cases. While assuming k! £ [0, 100] (higher redshifts) 
or fc' = 100 the best-fitting kinetic luminosity breaks shift 
significantly to higher values by ~ 1.5 decades. The other 
parameters do not seem to undergo any significant change. 
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Figure 6. Joint confidence intervals for the independent- z fits of Model S in z\ redshift range are shown [zi ^ 0.3). 68.3 per cent 
(solid, black), 95.4 per cent (dotted, red) and 99.7 per cent (dashed green) contours, based on AC statistics (see H4.3H . are shown. The 
best-fitting parameters are consistent with the data at the 90 per cent confidence level (p — value = 0.234). 
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Figure 7. Joint confidence intervals for the independent-z fits of Model S in Z2 redshift range are shown (0.3 < Z2 ^ 0.8). 68.3 per cent 
(solid, black), 95.4 per cent (dotted, red) and 99.7 per cent (dashed green) contours, based on AC statistics (see N4.3P ■ are shown. The 
best-fitting parameters are consistent with the data at the 90 per cent confidence level (p — value = 0.639). 
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Figure 8. Joint confidence intervals for the independent-z fits of Model S in 2:3 redshift range are shown (0.8 < £3 ^ 2.0). 68.3 per cent 
(solid, black), 95.4 per cent (dotted, red) and 99.7 per cent (dashed green) contours, based on AC statistics (see H4.3H . are shown. The 
best-fitting parameters are consistent with the data at the 90 per cent confidence level (p — value = 0.925). 




Figure 9. Overlaid 90 per cent confidence intervals based on AC statistics ( H4.30 for the three redshift ranges considered (zi drawn in 
black, Z2 in blue, 23 in red) of the independent-z fits of Model S. 
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Figure 10. Joint confidence intervals for the independent- z fits of Model W in Z\ redshift range are shown {z% ^ 0.3). 68.3 per cent 
(solid, black), 95.4 per cent (dotted, red) and 99.7 per cent (dashed green) contours, based on AC statistics (sec £|4,3l L arc shown. The 
best-fitting parameters are consistent with the data at the 90 per cent level (p — value = 0.251). 
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Figure 11. Joint confidence intervals for the independent-2 fits of Model W in 22 redshift range are shown (0.3 < 23 ^ 0.8). 68.3 per 
cent (solid, black), 95.4 per cent (dotted, red) and 99.7 per cent (dashed green) contours, based on AC statistics (see !|4.3t . are shown. 
The best-fitting parameters are consistent with the data at the 90 per cent level (p — value = 0.584). 
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Figure 14. From top right: (a) the joint confidence intervals of the best-fitting Qb and p m for population of sources with fixed maximum 
lifetime of t maXm = 8.5 X 10 6 yr is presented, (b) with t maXm = 1.2 X 10 7 yr, (c) with i maXm = 6.8 X 10 7 yr, (d) with t maXm = 1.3 X 10 8 
yr, (e) with t maXm = 5.4 X 10 s yr, and (f) with the maximum source age of t maXm = 1.1 X 10 9 yr. In each case t maXm is drawn from 
log-normal distribution with 0"t maXm = 0.05 (see 8)4.2.31 1. The independent- z fits of the Model S in the lowest redshift range (z± ^ 0.3) 
is presented. For each sub-plot the solid (black) line corresponds to 68.3 per cent of the best-fitting parameters for the particular choice 
of fixed tma X , the dotted (red) line corresponds to 95.4 per cent, and the dashed (green) contours contain 99.7 per cent of the best fits 
based on the AC statistics. It is clearly visible that the assumption on the maximum allowed age of the sources will strongly determine 
the density of the environment in which sources evolve (predominantly), as well as their jet powers (less strongly). 



Indeed, the addition of protons in the relativistic jets 
would require these particles, similarly to electrons, to be 
acc elerated to the relativistic speeds. As already shown by 
e.g. Bell (1978), relativistic protons store much more energy 
than electrons typically do, and this in turn will give rise to 
the radio source kinetic luminosity. The argument may be 
turned around: to maintain the expected kinetic luminosity, 
the observed radio lobe luminosity is much lower for pro- 
ton dominated jets than for lightweight, electron-positron 
jets (see also KDA). However, ma ny authors favour proton- 
electron jets. Ilstomin fc Soil <|201ll ). for instance, draw a hy- 
pothesis that some low-luminosity AGN have proton domi- 
nated jets, whi ch are suppo r ted b y extracting energy from 
black hole spin. ISikora et all (|2005t l report that jets of radio- 
loud quasars are also most likely heavy due to the positron- 
electron kinetic energy b eing too small to support en e rgetic s 
of radio lobes. However, Bickncll, Wagner & Groves (2001), 
extending the work of lCelotti fc Fabian (|l993l ). show that in 



the case when jets are proton dominated the following cri- 
terion occurs: jmin — 10 2 . This is not required by positron- 
electron jets (see £|5.4.4[> . Indeed, if one considers the net 
effect of the two assumptions tested, that is the jet parti- 
cle content and the minimum Lorentz factor, one will notice 
that their effects may cancel out. Setting 7 m i n — 10 2 or more 
one ensures that only initially highly relativistic particles 
are included in the jet material and these are later being re- 
accelerated in the jet front; it requires stronger momentum 



flux to affect the nonrelativitic protons than relativistic ones, 
hence the change in kinetic luminosity in these two cases. 

We are unable to distinguish between purely electron- 
positron jets (k' — 0) and those that contain modest num- 
bers of protons (k' £ [0, 10]). Very heavy FR II jets (k' ~ 
100) consisting of an electron-proton plasma seem to be un- 
likely, unless the minimum Lorentz factor is set to large val- 
ues ( 7min 10 2 - 10 4 ). 



6 DISCUSSION 

6.1 The lifetimes of radio galaxies 

Our estimated maximum lifetimes of radio galaxies may 
seem to be rather low, especially given that sources 
of the observed ages of order of few xlO 8 yr are 
known (e.g. ~ 1.8 X 10 8 yr of B0319-45 4 reported by 
ISaripalli. Subrahmanvan fc Hunstead1ll994 ~ 1.4 x 10 8 yr 



in th e case of B0313-683 reported by Schoenmakers et al.l 
1998). However, our results are in agreement with other 
estimates of observed ages and total lifetimes of, specif- 
ically, 3CRR radio galaxies evaluated by either nu- 
merical s imulations or multi-frequency radio observa- 
tions (e.g. Alex ander & Leahy 1987k_ 



Machalski et al.ll2007l : 



IWang fc Kaiserl |200S| ; lO'Dea et al] 1200911 . Our resu l ts are 
also consistent with work of Bird. Martini fc Kaiserl (|2008r ) 
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based on an independent low redshift sample, who deter- 
mined the maximum lifetime of the galaxies to be 1.5 x 
10 r yr. Further, the ages estimated by Machalski et al. 
(2004a) reach 1 — 2 x 10 s yr, but only for giant radio sources, 
and for most normal radio galaxies they analysed, i.e. the 
ones with linear sizes smaller than 1 Mpc, ages few— 10 times 
smaller were estimated. As a sanity check we find that for 
our hardwired maximum head advance speed of 0.4c a source 
of 1 Mpc in size may live for the maximum of 4 x 10 6 yr, 
while to increase the total lifetime to 1.6 x 10 s yr, the source 
head advance speed would have to be of order of 0.01c. Such 
a low head advance speed may be treated as a lower limit 
because the jet must be supersonic during the sources life- 
time to maintain the FR II morphology, and also classical 
double source s olde r than the inferred age are not observed. 
The model of lFalld (|l99lh predicts that the source advance 
speed changes, i.e. decreases, with the age of the source. This 
would suggest that sources that initially have high maximum 
head advance speed, slow down as they get older and hence 
they reach older ages and smaller w max are measured (e.g. 
iMachalski et al.ll2007l V 

The radio samples we are currently using contain the 
most powerful FR II radio galaxies and quasars, and the 
bias originating from the flux limit is particularly strong at 
high redshifts. This causes us to observe only the you ngest 
and most powerful sources. iBlundell fc Rawlingsl (|l999h sug- 
gested an explanation of the noticeable decrease with red- 
shift in source lifetimes. The term 'Youth-Redshift degen- 
eracy' has been coined to describe the effect. Due to the 
luminosity-redshift degeneracy originating from the radio 
flux limits, at high redshifts one can observe only the most 
powerful objects. They seem to be also of smaller linear 
sizes as compared to the radio galaxies in t he local Universe 
llKapahi. Kulkarni fc Subrahmanval Il987l ; iBarthel fc Milevl 
Il988h . The leading theoretical models of radio galaxy time 
evolution (e.g. KDA, BRW) predict a short initial phase 
of the radio source growth where the radio lobe luminos- 
ity peaks (up to a ~ few kpc), after which the radio lobe 
luminosity decreases as radio source ages and grows big- 
ger. Therefore, to be able to produce linear sizes and ra- 
dio lobe luminosities as observed at high redshifts, the 
sources must be younger than their low redshift counter- 
parts. Based on the 3CRR and BRL samples our results 
indeed suggest a trend of the total radio source lifetimes 
decreasing with redshift. The best-fitting parameters found 
for the combined-2 fits of the Models S and W indicate 
n t 6 [—3.5, —4.5] , slightly stronger a s compared to the re- 
sults obtalned_bvJ^^r^le --- eVaJ (| 1998h who found n t = —2.5 
and IWang fc Kaiser! (120081 ) who report n t = -2.4. How- 
ever again, due to the degeneracy between t maXm and p m 
the strength of the possible evolution is not strongly con- 
strained as it may be easily compensated by the redshift 
evolution of the central densities of the galaxy clusters or 
the kinetic luminosity break. Interestingly, a constant maxi- 
mum lifetime for all redshifts is allowed, i.e. n t ~ 0, but this 
would require the cosmological evolution of p m to be unfeasi- 
bly strong (see fl6.2p . Unfortunately, at the current stage we 
are unable to verify whether the decrease in the maximum 
lifetimes of these radio sources is genuine, or whether this is 
just an observational bias due to the flux limit. Deeper high 
redshift radio samples are needed to resolve this issue. 



6.2 Ambient gas densities 

The central gas densities estimated here are c onsistent with 
recent X-ray studies (e.g. iBelsole et al.ll2~007l: ICroston et aD 



l2008bl ; ISa"nderson, O'Sullivan fc Ponmanll2009l '). We observe 
that clearly higher densities are favoured for earlier epochs. 
Redshift evolution of the source environments is expected, 
at least, due to the expansion of the Universe. Invoking text- 
book physics it is known that 



H 2 (z) = ff?(fi M (l + zf + n k (i + zf + Oa), 



(16) 



where H is the Hubble parameter, Ho is the Hubble con- 
stant, and the parameters Qm, and S7k define the matter 
density, the vacuum density and the spatial curvature re- 
spectively. From the above equation one can deduce that 
the critical density evolves (p c ) with redshift as 



*>(*) = |% x [n M (i + z) 3 + n A ] 



(17) 



where G is the gravitational constant and assuming a flat 
Universe (fik = 0). Since the radio galaxies are often found 
in clusters and galaxy groups, one may expect that their 
immediate environments will undergo similar evolution to 
the one pred icted by models of cluster formation (see e.g. 
lArnaudll2005l . for a review). Structure formation models pre- 
dict that the mean dark matter density (the predominant 
constituent of galaxy clusters) scales as the critical den- 
sity of the Universe, that is again p (z) oc p c (z) oc H 2 (z). 
The critical density will, therefore, change between redshifts 
z n 2 and z n i by a factor of p c (z n 2) / p c (z n i) ■ If z n i = and 
z n 2 = 2 then the density would be expected to be higher 
at z n 2 by a factor of z n 2/z n i = 8.8, equivalent to ~ 0.94 in 
logarithmic scale, and to the redshift dependence as approx- 
imately (1 + z) 2 . Simil ar trends have been recently found by 
IPoggianti et~al1 (|2010l ). and as they point out it ca n be also 



deduced from previous o bservational studies fe.g. iDresslei 
1980t iDressler et al. I ll997l) . On the other hand. lO'Dea et al 



l|2009h do not observe any density evolution up to z = 0.5. 
Our best-fits for the combined-z fits of the tested models 
(Table 3} clearly show that n r is degenerate and strongly 
dependent on nt and n q ; the value expected by the galaxy 
cluster cosmological evolution models and the effect of the 
expanding Universe is possible for certain values of n t and 
n q , but is not unique. Some may argue that environments 
of FR II sources may depend on the cosmological epoch. 
That is, beside the argument of the evolution following the 
expansion of the Universe, at low redshifts weak clusters 
and galaxy groups are favoured by FR II sour ces, while rich 
galaxy clusters are preferred at high redshifts (|Deltorn et al.l 
Il997l ; IZirbellll997l ; IWold et al.ll2000D . If this is the case the 
effect will be convolved in the n r parameter resulting in evo- 
lution stronger than the expected n r ~ 2.0 — 2.5. The possi- 
ble values of n r and their consequences are discussed further 
in 3Q1 

As has already been shown in H5.4l the spread of the log- 
normal distribution of the density has a rather small effect 
on the results. This is highlighted by the fact that when a 
delta function is employed the final results are still in good 
agreement with the observed samples. We conclude that at 
the current stage we are unable to constrain the shape of the 
population ambient density distribution. Further, the mean 
central density value will depend on the assumed a and ft, 
and hence if our assumptions for these values are incorrect 
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Figure 15. 68.3 per cent (solid, black), 95.4 per cent (dotted, red) and 99.7 per cent (dash-dotted, green) joint confidence intervals 
(based on AC statistics) of all the searched parameters of the combined- z fits of the Model S. Parameters nt,n q and n T quantify the 
strength of the redshift evolution of t m ax m , Qb and p m respectively. 
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Figure 16. 68.3 per cent (solid, black), 95.4 per cent (dotted, red) and 99.7 per cent (dash-dotted, green) joint confidence intervals 
(based on AC statistics) of all the searched parameters of the combined- z fits of the Model W. Parameters nt,n q and n r quantify the 
strength of the redshift evolution of t m ax m , Qb and p m respectively. 
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Figure 17. Kinetic luminosity functions generated with the best-fitting parameter sets of the independent-z fits of Model S (solid 
line, green) and Model W (solid line, black) for all redshift ranges, and combined- z fits of Model S (blue) and Model W (orange). The 
combined- z fits of Model S and Model W show their best fits (solid blue, solid orange respectively), and two restricted cases where n T 
parameter is set to n T = 2.5 (dashed) and n r = 4.0 (dash-dot-dot-dotted). 



one may simply scale the fitted densities with the (pofflj) 
relation. We currently assume that there is no redshift evo- 
lution of the core radius or the density profile slope. If there 
is any intrinsic change this will be subsumed by the mean 
density value, possibly resulting in a stronger trend with 
redshift. 



6.3 Kinetic luminosities of powerful FR II sources 

The kinetic luminosity functions (KLFs) of the radio source 
populations generated with the best-fitting parameters for 
all redshift ranges and models are presented in Figure 1171 
while in Figure [18] the kinetic luminosity distributions and 
their corresponding underlying parent populations from 
which the observed samples originate are plotted for each 
redshift range, size bin and model. There are two main con- 
clusions. Firstly, the higher the redshift the more power- 
ful the jets are that are included; at z < 0.3 the high- 
est Q are of order 10 40 W (Model S), while in the high- 
est redshift range, Z3, the most powerful sources seem to 
have jet luminosities reaching 10 42 W (Figure [17)) . Secondly, 
the larger the radio sources are the more powerful AGN 
on average they host (Figure [18)) . The second effect seems 
to be much milder than the dependence on redshift, but it 
is present; a similar pattern has been a lready observed by 
iGopal-Krishna. Wiita fc Saripallil l) 19891 ). 



'erences between results of Model S and W 



6.3.1 

As can be seen from Figure [17] and Figure 1181 although 
the parent populations of Model S and W are different, the 
predicted flux-limited samples of the two models are very 
similar. More interestingly, one would expect kinetic lumi- 
nosities of sources of FR II type to be at least 10 37 W. In- 
deed, although we simulated sources of kinetic luminosities 
of 10 36 W, these do not contribute to the observed popula- 
tions. Some consideration is also needed of the lowest red- 
shift bin best fits for Model W, where the kinetic luminosity 
break reaches somewhat low values and the slope is flatter 
compared to the other redshift ranges. The kinetic lumi- 
nosity distribution function of Model W has been used by 
Willott et alJ (|200lh as their initial kinetic luminosity dis- 
tribution function but for the most powerful sources only. 



In their work the population of radio galaxies and radio- 
loud quasars was divided based on the emission line char- 
acteristics rather, that is high- and low-excitation sources 
(HEG and LEG respectively), than the historical FR I and 
FR II classes. HEG include powerful FR II and the most 
powerful FR I sources, while LEG include most FR Is and 
low-luminosity FR II objects. As mentioned before Model 
W may not be able to reproduce th e faint end of the dis- 
tribution. Contrary to IWillott et alJ ()200lh we do include 
low luminosity FR lis and perhaps this is the cause of the 
observed shift of Qb and a s . 



6.3.2 Bending power-law as the initial distribution of 
kinetic luminosities 

The thin, dotted lines in Figure [18] represent the possible 
parent populations generated according to the assumed dis- 
tribution functions. The flux limit of the radio catalogues 
used is very high, and hence only the most powerful sources 
are observed; as estimated by our simulations their kinetic 
luminosities do not drop below 10 37 — 10 38 W. In such a 
strong flux limit regime either of the kinetic luminosity dis- 
tribution functions (Models S and W) will shape the faint 
power end of the luminosity functions on the basis of a good 
fit to the high power end. However, with the lack of observed 
sources with Q < 10 38 W the faint end of the function is 
highly unconstrained especially at higher redshifts where, 
given currently used radio samples, only the most power- 
ful sources are considered. To investigate this issue we at- 
tempted to fit pre- and post- kinetic luminosity break slopes 
separately. To do this, we have tested the so-called curved 
power-law distribution 



i>{Q) dQ = V* 



Q 



+ 



dQ. 



(18) 



where a s \ describes the pre-break slope of the func- 
tion, and « S 2 the post-break slope. To avoid introduc- 
ing more free parameters we decided to set the high 
power end of the kinetic luminosity distribution function 
to q B 2 = 2.2 (after the radiative lu minosity functions of 
Hopk ins. Richards k, Hernauisdl2007h . and allowed the pre- 
break slope to vary. The results, while assuming the curved 
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power-law as the initial distribution function of source ki- 
netic luminosities, are in good agreement with those found 
with the Model S. The faint end of the kinetic luminosity 
function is well constrained in the lowest redshift range, but 
since now it is not dependent on the high power end of the 
function it becomes highly unconstrained at higher redshifts. 
It is necessary, therefore, to include lower luminosity radio 
samples in our analysis to investigate the faint end of the 
kinetic luminosity functions. 



6.3.3 Cosmic downsizing 

The concept of cosmic downsizing, as an anti-correlation 
between the stellar mass of galaxie s and the e p och o f 
their formation, was introduced by ICowie et all l|l996h . 
Since then, the term has been broadened to include the 
AGN evolution that has manifested similar anitihierarchi- 
cal trends. The suggestion that the low-luminosity AGN 
X-ray luminosity functions seem to pe ak at lower red- 
shifts than those of the luminous q uasars (|Ueda et alj|2003l ; 
lHasinger. Mivaii fe Schmidd l2005). has been interpreted as 
implying black hole anitihierarchical growth, and hence 
seems to b e analogous to the st ar-forming galaxy down- 
sizing (e.g. IMerloni fc Heinzl 120081 ). Furthermore, since the 
growt h of the black hole and its ki netic luminosity may be 
linked. IMerloni fc Heinz! (|2008l ) and lCattaneo fc Best! <|2009l ) 
suggested that the kinetic energy density of these active 
galaxies undergoes a similar cosmological evolution such 
that the most p o werfu l jets are found at earlier epochs. 
IMerloni fc Heina l)2008l ) based their findings on the anal- 
ysis of radio core kin ematic luminosity functions, while 
ICattaneo fc Bestl l|2009t) analysed imprints of the extended 
radio emission on the X-ray emitting plasma. It should not 
be surprising, therefore, if such an effect is seen in our re- 
sults. The results of the independent- z fits already suggest 
a redshift evolution of the kinetic luminosity break. More- 
over, the results of the combined- z fits of the fitted models 
confirm that this trend with redshift may be rather strong. 
Of course, the strength of this cosmological evolution de- 
pends on the assumed evolution of the sources' maximum 
lifetime and/or environment in which sources grow due to 
previously discussed degeneracies. One of the most plausible 
cases, where the galaxy cluster central environments evolve 
with the expanding Universe and according to the galaxy 
cluster evolutionary models, would require n r = 2.0 — 2.5, 



r](a)mc 



(20) 



and this will lead to n q = 9.5 and n t 



-3.0. Given that 



Qb (z = 0) = 2 x 10 W here, such a strong redshift evolu- 
tion would imply Q B (z = 2.0) ^ 6.8 x 10 42 W. This value is 
rather high, the highest kinetic luminosities of FR II radio 
sources estimated in other works are of order of 10 41 W. 

To give an estimation on the black hole mass which 
would produce such a kinetic luminosity one needs to know 
the black hole r adiative, e(a), and jet, 77(a), efficienc ies. 
Specifically, after IMartmez-Sansigre fc Rawlingj (|201lh . it 
is assumed here that the bolometric luminosity of an AGN 
accreting at the rate m is 



Lboi = e(a)ihc 



(19) 



where c is the speed of light. Furthermore, the AGN's kinetic 
luminosity depends also on the accretion rate rh and is given 
by 



Since the Eddington ratio is defined as AEdd = iboi/iEdd 
one can easily derive the expression for black hole mass, 
which is 



Mb 



0.8 x 10" 



e(a) 



?7(a)A E 



-Q 



(21) 



in the units of Mg. For non-spinning black holes, that is 
assuming that these radio sources power their jets via ac- 
cretion at their Eddington limits (AEdd ~ 1-0), and so 
the spin a = 0, the most powerful objects inferred here 
at z ~ 2.0 (Q B = 6.8 x 10 42 W) would host black holes 
of masses of ~ 2.7 x 10 14 Mq. For this case, the jet effi- 
ciency is found to be within 77(a) € (10 _5 ,10 -3 ) depend- 
ing on the theoretical model assumptions, and the radia- 
tive efficiency is e(a) ~ 0.05 (see references compiled by 
IMartmez-Sansigre fc Rawlingsll201ll ). Similar results are ob- 
tained for black holes of low spins (e.g. for a — 0.4 the ef- 
ficiencies are 77(a) « 0.001 and e(a) f» 0.05 which imply 
~ 2.7 x 10 13 M Q ). However, assuming that the black holes 
are of high spins (analogically setting 77(a) « 0.5 and e(a) « 
0.1 which imply black hole spin of a > 0.95) the black holes' 
masses would need to be ~ 1.0 x 10 11 Mq in order to produce 
quoted Qb ■ The most massive AG Ns' black holes are found 
to not exceed ~ 10 10 M q (e.g. IkcLure fc Dunlod [2004 



iNataraian fc Treistedl2009l ; ISteinhardt fc Elvisll2010l . l201ll ). 
Yet, equivalently to the independent-z fits, at these high 
redshifts the Qb values are of such small number density 
that the highest kinetic luminosities that contribute to the 
population do not reach 10 42 W (Figure [TTJl . In such a case, 
assuming maximum kinetic luminosity of Q = 1.0 x 10 42 W 
the non-spinnig black hole as denoted earlier would need to 
have a mass of ~ 3.8 x 10 13 Mq, and the high spin black 
hole would need to have a mass of ~ 1.5 x 10 10 M Q . 

Energy extraction from black hole spin is an attrac- 
tive way to explain high jet luminosities and efficien- 
cies (e.g. iPenrose fc Flwdlll97ll ; iBlandford fc ZnaiekHl977l ; 
iRees et al.lll982l ; lMeierlll999h . After iMeierl (120021 ) we note 
that Q oc j 2 B 2 ol (MBH/MQ) 2 , where j denotes the black hole 
spin, and B po \ is the poloidal magnetic field. In such a case, 
for a given kinetic luminosity the higher the spin is, the 
less massive the black hole must be. However, one must 
bear in mind that the strength of the poloidal magnetic 
field will have a direct effect on the efficiency of the black 
hole rot ational ener gy extraction (for discussion on this is- 
sue see lMeierll2002l ). A combination of accretion and spin 
may play a significant role in producing the kinetic lumi- 
nosities estimated here. However, our results suggest that 
also environments and lifetimes of radio galaxies have a di- 
rect influence on the radio luminosity of large scale struc- 
tures of radio galaxies. It is worth considering, therefore, 
whether it is possible to reduce the redshift evolution in Qb 
based on the interplay between the environments and radio 
galaxy lifetimes in our results. For instance, the Qb cosmo- 
logical evolution of n q = 7.0 will yield milder evolution of 
the radio source maximum lifetimes, n t = —1.0, and the 
trend with redshift of the galaxy cluster central densities 
can be as strong as n r = 10.5 (Model S) or even n T — 11.5 
(Model W). With this pattern of cosmological evolution 
the kinetic luminosity break would reach ~ 1 x 10 41 W at 
z — 2.0, which further implies maximum black hole masses 
of 10 10 Mq if the sources accrete close to their Eddington 
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Figure 18. Functions showing kinetic luminosity distributions (normalized for each redshift range separately) of the observed, that 
is flux limited, sources (solid lines) and the corresponding underlying parent populations from which the observed samples originate 
(dotted lines). In each plot, for each redshift range and model considered, kinetic luminosity distributions for separate linear size bins 
(as used in fitting) are plotted. The exact radio galaxy linear sizes (D) falling into each size bin are detailed in corresponding panels. 
The distributions change for different redshifts as well as for different radio galaxy linear sizes, where larger radio sources on average 
seem to be also more powerful. Note that these distributions of parent populations depend strongly on the slope a s which for z > 0.3 
is associated large uncertainties. Here, for pictorial purposes, only the best-fits for each redshift range are presented. Colours denoting 
linear sizes of the radio galaxies are labeled. 
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luminosity, and jets are predominantly accretion powered. 
However, this causes n r to be unfeasibly strong. One may 
get more realistic results while assuming n r = 4.0, as then 



The duty cycle, <5(Mbh, z), is defined as 



our results show that n q = 8.0 and n t 



-2.5, which leads to 



5(M BH ,z) = 



QbO = 2.0) 6.6 x 10 41 W, t m ax m (z = 2.0) 3.9 x 10 6 yr, 
and Pin (z = 2.0) 4 x 10 22 kgm- 3 (Model S). The red- 
shift evolution of the central densities may seem to be rather 
high, since at n r = 4.0 it would exceed the effect from the 
Universe's expansion; however, as discussed in H6.2I such a 
strong evolution may be realistic if some intrinsic change in 
the environments, or other effect such as evolution of core 
radius, take place. Alternatively, radio sources may be al- 
lowed to reach older ages than the ones found with the best- 
fitting parameters (but kept within the confidence intervals 
as found in this work). To lower the kinetic luminosities the 
sources would have to be younger or reside in higher den- 
sity environments than our best-fitting parameters suggest, 
or undergo a combination of the two effects. For instance, if 
one attemp ts to match the kinetic lumin osities and ages es- 
timated bv lRawlings fc Saundersl (|l99lh for 3CRR sources, 
the lifetimes would have to be allowed to reach 10 s yr, and 
hence sources would have to reside in very high density en- 
vironments. Finally, there is a possibility that, contrary to 
what we assumed so far, the intrinsic source parameters are 
not independent. Scenarios where Q and t, and/or t and Rt 
are correlated are sometimes considered, but are beyond the 
scope of our current work. 

It is necessary to point out that as the cosmological evo- 
lution of the source parameters is certainly real and is not 
an effect of an observational bias due to high flux limit of 
the observed samples, the location of the kinetic luminosity 
break can be only constrained in the highest redshift bin. 
From our results it seems that in the lower redshift bins the 
Qb must be at the lower end of the luminosities that still 
contribute to reconstruction of the observed samples, but its 
exact location is uncertain (Figure [TSJ .In addition, one may 
notice in Figure [18] that the independent-z fits predict lower 
number density of sources than the combined- z fits at inter- 
mediate and high redshifts. The reason for this behaviour 
lies in tying the a s exponent of Eqn.[8]and Eqn.[9]in the case 
of the combined-z fits, where a s is required to stay the same 
for all redshifts (unlike for the independent-z fits where the 
exponent may vary freely between redshift ranges) . The best 
fits of both model S and W of the combined-z fits predict 
much flatter slopes of the kinetic luminosity distributions 
than the independent-z fits. This will cause the existence of 
low radio luminosity sources that are not observed due to 
the high flux limit of the samples, but are unveiled in the 
KLFs. Both of these issues this may be verified with larger 
samples with lower flux limits than are currently being used 
in this study. 

6.4 Constraints on AGN duty cycles 

The black hole masses of massive elliptical galaxies, 
which are ty pical hosts of powerful r adio galaxies, are 
known Je.g. IWolfe fc Burbidgel 1 197(1 IWu. Liu fe Zhangj 



120021 ; iDunlop et all 120031 ; [Gadotti fc Kaufmannl 12009ft . 
Hence, once number densities and kinetic luminosities of the 
sources have been estimated, one can set constraints on the 
AGN duty cycle which defines the fraction of black holes 
being active at a particular cosmological epoch. 



0(MactivoBH, z) 



BH, Z) 



(22) 



where </>(M ac tiveBH, z) is the space density of active AGN 
with mass M and at the redshift z, while <^(Mbh,z) anal- 
ogously defines the space density of all AGN (see also 
iKellv et al.ll2010ft . Since in our work other types of radio 
loud AGN are excluded, e.g. FR Is and blazars, we shall 
refer to the duty cycle of FR II radio galaxies only in this 
discussion; this will be denoted by 5frii(A/bh, z). 

A typical FR II source can be assumed to have 
a black hole mass of ~ 10 9 M©. If we assume that 
these are the most typical sources during each epoch 
(i.e. there is no significant mass evolution), then the typ- 
ical kinetic luminosity of these black holes is around 
the value of Qb of the considered cosmological epoch. 
For the lowest redshifts, with Q B =± 10 39 W, the 
space number density is approximately 0(M act ivoBH, zi) ~ 
10" 7 dN Mp C - 3 Alog 10 (Q)- 1 (FigureE}; 0(M act iveBH, z) is 
found analogously for the other two redshift bins. These val- 
ues are com pared with the local black hole mass function s 
compiled by IShankar. Weinberg fc Miralda-Escudel {2009), 
which give approximately 4>(Mbb = 10 9 Mq,z = 0) = 
10" 4 dN Mpc" 3 Alog 10 (Af B H) _1 . Assuming that, around 
10 9 M , the average fractional accretion rate and efficiency 
of conversion to kinetic power are independent of black 
hole mass, we can infer that Alog 10 (Q) = Alog 10 (AfBH). 
Then, neglecting redshift-evolution of the black hole mass 
function, we find duty cycles of <5frii(A4bh, zi) ~ 10 -3 , 
5frii(M BH ,Z2) ~ 10" 4 - 5 , and 5 F rii(Mbh, z 3 ) ~ 10" 5 5 . Al- 
though in our analysis we use ranges of redshifts, the mea- 
surement is instantaneous, and hence can be compared to 
other results for a single redshift; here our measurement will 
be an average for the redshifts within each range we consider. 

We can compare the estimated duty cycles with the 
maximum lifetimes of radio galaxies obtained from our fits 
to determine whether all or just a fraction of 10 9 Mq black 
holes are likely to spend some time as a powerful FR II ra- 
dio galaxy, where we define 'powerful' to mean at the break 
luminosity, i.e. the most powerful 'typical' objects in the ob- 
served epochs. From our fits, we found that the sources live 
for a maximum of approximately 10 7 ' 8 yr in our lowest red- 
shift range, and since Azi = 10 9 ' 5 yr, then each radio galaxy 
is on for 1.6 per cent of this epoch, compared to a duty cy- 
cle of 0.1 per cent. If we assume that an AGN goes though 
the FR II stage only once, if at all, during its lifetime, then 
around 10 per cent of 'typical' radio galaxy hosts are likely 
to spend time as a powerful FR II in the most recent epoch. 
Moreover, if we assume that these galaxies can go through 
their FR II episodes multiple times throughout this epoch, 
then the fraction will be even smaller. For the higher red- 
shifts the inferred duty cycles are smaller as compared to 
the local Universe. From the fits, the radio sources at higher 
redshifts are found to live for a shorter time; the lifetimes of 
10 72 - 10 7 5 yr for redshift range z 2 and 10 B 6 - 10 72 yr for 
Z3 are found, while the time spanned by each epoch remains 
close to Azi. However, the reduction in lifetimes is not suf- 
ficient to explain the drop in duty cycle at high redshifts: 
the most powerful FR lis become rarer, and a smaller frac- 
tion of potential FR II hosts spend time as powerful FR lis 
(perhaps < 1 per cent for z%, and < 0.1 per cent for Z3). 
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These apparent reductions in duty cycle may result 
from our assumption of no evolution in the number den- 
sity of 10 9 M0 black holes from the present epoch. In 
fact, based on th e evolu tion of X-ray luminosity functions, 
Merloni fc Heinzl (|200dh (see also iNataraian fc Treisteil 



modified models of lBarai fc Wiital (|2006l , 120071 ), might yield 
somewhat different results. 
We find that: 



2009) estimate that the black hole number density of SMBHs 



10" 



to 



declines from (j){M m = 10 9 M Q ,z = 0.3) 
<f>(M m = 10 9 M Q , 2 = 2.0) ~ 10" 5 - 8 . This reduction in num- 
ber density at higher redshifts could account for most of the 
apparent drop in duty cycle which we infer from our fits. 
It is also important to note that since the break luminos- 
ity of the kinetic luminosity function increases by 2 decades 
from z\ to 23, the active fractions quoted here are tracing 
significantly different populations in terms of kinetic power, 
and if lower-power radio galaxies are included (e.g. the FR I 
population) the total duty cycle of radio-loud AGN may be 
much greater than the values given here. 

It is interesting to compare our estimates of powerful 
radio galaxy duty cycles with estimates of AGN duty cy- 
cles based on radiative luminosity functions. In the higher 
redshift ranges, our duty cycle estimates are at least a 
decade lower than those inferred for 1 9 Mq black holes by 
IShankar. Weinberg fe Miralda-Escudl (120091 ) and the lower 
bounds estimated by iKellv et al. I (|2010l)~ This discrepancy 
may be understood in terms of the radio loud fraction of 
AGN, since for the high kinetic powers observed at higher 
redshifts, the accretion rates are expected to be large frac- 
tions of the Eddington rate and hence correspond to the 
underlying quasar population, which is known to be pre- 
dominantly radio-quiet. Hence at these redshifts, the lower 
duty cycle inferred from kinetic luminosities may reflect 
the small fraction of high accretion-rate, radio-loud objects. 
However, at lower redshifts, the situation turns around, 
and our inferred duty cycle for 10 9 Mq may be up to a 
decade larger than that estimated from radiat i ve lum inosi- 
ties by IShankar. Weinberg fc Miralda-Escudel (|2009T l. This 
difference may reflect the evolution of the FR II kinetic 
break luminosity to lower values, which may correspond to 
a lower typical black hole mass at low redshifts than the 
10 9 M Q we assume here, for which the du t y cycl es predicted 
by IShankar. Weinberg fc Miralda-Escudel |2009) are higher. 
An alternative possibility is that the reduction in break lu- 
minosity corresponds to a change to a lower accretion rate 
(as a fraction of the Eddington rate), where a kinetically- 
dominated accretion mo de dominates over radi atively effi- 
cient accretion (e.g. see iMerloni fc Heinzl 120081 . and refer- 
ences therein). 



7 SUMMARY 

We perform multidimensional Monte Carlo simulations to 
investigate fundamental source parameters (such as kinetic 
luminosity, age and ambient density), and their possible 
trends with redshift of powerful radio galaxies and radio- 
loud quasars of FR II morphology. We present a new 
method, based on radio luminosity functions, for estimating 
the parent populations of complete, flux limited radio sam- 
ples. In our analysis we use the KA/KDA semi-analytical 
model of FR II type radio source time evolution, and we 
stress that other models, such as those of BRW, MK, or the 



(ih) 



(v) 



(vi) 



(vii) 



(viii) 



The total lifetimes of the radio galaxies are found to 
be few x 10 7 yr at low redshift and decrease for earlier 
epochs. This is in agreement with independent studies 
on the 3CRR radio sources, but may be specific to 
these most powerful radio galaxies and quasars. With 
the current sample and its strong flux limit we are 
unable to draw final statements about the evolution of 
the lifetimes of the sources which is suggested by our 
current results. 

Our results suggest cosmological evolution of one or 
more source parameters. In particular the mean den- 
sity of the immediate source environments (or the p af 
parameter if we allow core radius or j3 to change with 
redshift) is found to undergo evolution with redshift; 
the hypothesis that there is no evolution is ruled out 
with probability of 99 per cent. 

The central density of the FR II environments is 
found to undergo redshift evolution of approximately 
(1 + z) 4 ; evolution stronger than the one expected from 
the Universe expansion is possible if we consider addi- 
tional effects from change of environment or core radius 
evolution. 

The function describing the initial distribution of ki- 
netic luminosities modelled by the Schechter function 
(or its modification) or smoothly broken power-law is 
favoured to simple power-law; however, the hypothesis 
of an unbroken power-law distribution cannot be ruled 
out at a confidence level of more than 95 per cent. 

The estimated kinetic luminosities are within 10 38 — 
10 41 W for FR II type sources, which is consistent with 
previous studies. We observe that the kinetic luminosi- 
ties depend on the cosmological epoch and the linear 
size of the source, where larger in linear size as well as 
higher redshift sources are more powerful. 

The FR lis' kinetic luminosity function undergoes 
cosmological evolution of the break luminosities of at 
least (1 + z) 3 and may be as strong as (1 + z) 10 . The 
uncertainty originates from the strong degeneracy be- 
tween Qb, imax and p . Evolution stronger than (1+z) 9 
is rather unlikely since as the consequence the black 
hole masses of the most luminous FR II sources would 
have to be of > 10 11 Mq assuming that there is no 
strong spin powering of the jets. 

Our results suggest that, at least at high redshifts, 
FR II sources most probably accrete at moderate/high 
Eddington ratios and the black hole spin may play a 
significant role in the jet production, as both effects 
seem to be necessary to explain the high estimated 
kinetic luminosities of FR II sources at higher z. 

We estimated the duty cyles of FR II radio galaxies 
at the break in the kinetic luminosity function, find- 
ing them to decrease with redshift from ~ 10 -3 for 
zi < 0.3, to ~ 10" 4 - 5 at 0.3 < z 3 sC 0.8, to ~ 10" 5 5 at 
0.8 < Z3 2.0. The decrease in duty cycle at higher 
redshifts may be explained by a combination of the 
reduction in the lifetime of FR II radio galaxies to- 
gether with evolution in the number density of mas- 
sive black holes. The shift in kinetic luminosity break 
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to higher values also indicates an intrinsic change in 
the population of kinetic powers. Interestingly, at low 
redshifts the duty cycle of powerful FR lis exceeds that 
estimated for 10 9 Mq AGN based on radiative lumi- 
nosities. This difference can be explained if the typical 
black hole mass of FR lis shifts to lower masses at 
low 2. Alternatively, the low-z FR II population may 
become dominated by a kinetically dominated, radia- 
tively inefficient mode of accretion. 



ACKNOWLEDGMENTS 

ADK thanks Peter A. Curran for many valuable discussions. 
The authors thank Judith H. Croston, Thomas J. Maccarone 
and Andreas Papadopoulos for their useful and insightful 
comments on the manuscript. The authors thank the referee 
Paul Wiita for comments which improved the content of the 
paper. ADK acknowledges the Leverhulme Trust for the fi- 
nancial support. PU acknowledges funding from an STFC 
Advanced Fellowship. PU and ADK acknowledge funding 
from the European Community's Seventh Framework Pro- 
gramme (FP7/2007-2013) under grant agreement number 
ITN 215212 "Black Hole Universe". This work made an ex- 
tensive use of the Iridis Compute Cluster maintained by 
the University of Southampton, Southampton, UK, and, in 
the later stages of the project, of the Sciama High Perfor- 
mance Compute Cluster maintained by the University of 
Portsmouth and SEPNet (South-Eastern Physics Network) 
UK. 



REFERENCES 

Alexander P., Leahy J.P., 1987, MNRAS, 225, 1 

Allington-Smith J.R., Ellis R., Zirbel E.L., Oemler A.Jr., 
1993, ApJ, 404, 521 

Alshino A., Ponman T., Pacaud F., Pierre M., 2010, MN- 
RAS, 407, 2543 

Arnaud M., 2005, 'Background Microwave Radiation and 
Intracluster Cosmology', eds. Melchiorri F. & Rephaeli Y., 
Italian Physical Society Series, p. 77 

Ashkarian T.G., Longair M.S., 2000, MNRAS, 311, 846 

Barai, P., Wiita P.J., 2006, MNRAS, 372, 381 

Barai, P., Wiita P.J., 2007, ApJ, 658, 217 

Barthel, P.D., 1989, ApJ, 336, 606 

Barthel, P.D., Miley G.K., 1988, Nature, 333, 319 

Baum S.A., Heckman T.M., van Breugel W., 1992, ApJ, 
389, 208 

Baum S.A., Zirbel E.L., O'Dea CP., 1995, ApJ, 451, 88 
Bell A.R., 1978, MNRAS, 182, 443 

Belsole E., Worrall D.M., Hardcastle M.J., Croston J.H., 
2007, MNRAS, 381, 1109 

Best P.N., Bailer D.M., Longrair M.S., Riley J.M., 1995, 
MNRAS, 275, 1171 

Best P.N., Peacock J.A., Brookes M.H., Dowsett R.E., 
Rottgering H.J.A., Dunlop J.S., Lehnert M.D., 2003, MN- 
RAS, 346 1021 

Best P.N., Rottgering H.J.A., Lehnert M.D., [BRL], 1999, 
MNRAS, 310, 223 

Best P.N., Rottgering H.J.A., Lehnert M.D., 2000, MN- 
RAS, 315, 21 



Bicknell G. V., Wagner S.J, Groves B., 2001, ASPC, 250, 
80 

Bird J., Martini P., Kaiser C, 2010, ApJ, 676, 147 
Blandford R.D., Rees M.J., 1974, MNRAS, 169, 395 
Blandford R.D., Znajek R.L., 1977, MNRAS, 179, 433 
Blundell K.M., Rawlings S., 1999, Nature, 399, 330 
Blundell K.M., Rawlings S., Willott C.J., [BRW], 1999, 
ApJ, 117, 677 

Blundell K.M., Fabian A.C., Crawford C.S., Erlund M.C., 

Celotti A., 2006, ApJ, 644, 13 
Bohringer H., Voges W., Fabian A.C., Edge A.C., Neumann 

D.M., 1993, MNRAS, 264, 25 
Brent R.P., 1973, 'Algorithms for Minimization without 

Derivatives', Prentice-Hall, Inc., Englewood Cliffs N.J. 
Carilli C.L., Owen F.N., Harris D.E., 1994, AJ, 107, 480 
Cash W., 1976, A&A, 52, 307 
Cash W., 1979, ApJ, 228, 939 
Cattaneo A., Best P.N., 2009, MNRAS, 395, 518 
Celotti A., Fabian A.C., 1993, MNRAS, 264, 228 
Cowie L.L., Songalia A., Hu E.M., Cohen J.G., 1996, AJ, 

112, 839 

Croton D.J., Springel V., White S.D.M., De Lucia O, 
Frenk C.S., Gao L., Jenkins A., Kauffmann G., Navarro 
J.F., Yoshida N., 2006, MNRAS, 365, 11 

Croston J.H., Hardcastle M.J., Harris D.E., Belsole E., 
Birkinshaw M., Worrall D.M., 2005, ApJ, 626, 733 

Croston J.H., Hardcastle M.J., Birkinshaw M., Worrall 

D. M., Laing R.A., 2008a, MNRAS, 386, 1709 

Croston J.H., Pratt G.W., Bohringer H., Pointecouteau 

E. , Ponman T.J., Sanderson A.J.R., Temple R.F., Bower 
R.G, Donahue M., 2008b, A&A, 487, 431 

Curran P.A., Evans P.A., de Pasquale M., Page M.J., van 
der Horst A.J., 2010, ApJL, 716, 135 
Daly R.A., 1995, ApJ, 454, 580 
Daly R.A., 2011, MNRAS, 414, 1253 

Deltorn J.-M., Le Fevre O., Crampton D., Dickinson M., 

1997, ApJ, 483, 21 
Dressier A., 1980, ApJ, 236, 351 

Dressier A., Oemler A.Jr., Couch W.J., Smail I., Ellis R.S., 

Barger A., Butcher H., Poggianti B.M., Sharpies R.M., 

1997, ApJ, 490, 577 
Dunlop J.S., Peacock J.A., 1990, MNRAS, 247, 19 
Dunlop J.S., McLure R.J., KukulaM.J., Baum S.A., O'Dea 

CP., HughesD.H., 2003, MNRAS, 340, 1095 
Dunn R.J.H., Fabian A.C., 2004, MNRAS, 355, 862 
Eales S.A., 1985, MNRAS, 217, 149 
Falle S.A.E.G., 1991, MNRAS, 250, 581 
Fanaroff B.L., Riley J.M., 1974, MNRAS, 167, 31 
Gadotti D.A., Kaufmann G., 2009, MNRAS, 399, 621 
Gopal-Krishna, Wiita P.J., 2000, A&A, 363, 507 
Gopal-Krishna, Wiita P.J., 2001, A&A, 373, 100 
Gopal-Krishna, Wiita P.J., Saripalli L., 1989, MNRAS, 

239, 173 

Grimes J.A., Rawlings S., Willott C.J., 2004, MNRAS, 349, 
503 

Hardcastle M.J., Croston J.H., 2010, MNRAS, 404, 2018 
Hasinger G., Miyaji T., Schmidt M., 2005, A&A, 441, 417 
Hill G.J., Lily S.J., 1991, ApJ,367, 1 

Hopkins P.F., Richards G.T., Hernquist L., 2007, ApJ, 654, 
731 

Istomin Ya.N., Sol H., 2011, A&A, 527, 22 
Jarvis M.J., Rawlings S., MNRAS, 2000, 319, 121 



30 A.D. Kapniska, P. Uttley & C.R. Kaiser 



Jarvis M.J., Rawlings S., Willott C.J., Blundcll K.M., Ealcs 

S., Lacy M., 2001, MNRAS, 327, 907 
Jarvis M.J., Teimourian H., Simpson C, Smith D.J.B., 

Rawlings S., Bonfield D., 2009, MNRAS, 398, 83 
Kaiser C.R., Alexander P., [KA], 1997, MNRAS, 286, 215 
Kaiser C.R., Alexander P., 1999, MNRAS, 305, 707 
Kaiser C.R., Best P.N., 2007, MNRAS, 381, 1548 
Kaiser C.R., Best P.N., 2008, MNRAS, 384, 1742 
Kaiser C.R., Dennett-Thorpe J., Alexander P., [KDA], 

1997, MNRAS, 292, 723 
Kapahi V.K., 1989, ApJ, 97, 1 

Kapahi V.K., Kulkarni V.K., Subrahmanya C.R., 1987, 
JApA, 8, 33 

Kawakatu N., Kino M., Nagai H., 2009, ApJ, 697, 173 
Kelly B.C., Vestergaard M., Fan X., Hopkins P., Hernquist 

L., Siemiginowska A., 2010, ApJ, 719, 1315 
King I.R., 1972, ApJL, 174, 123 
Kino M., Takahara F., 2004, MNRAS, 349, 336 
Kirk J.G., Guthmann A.W., Gallant Y.A., Achtenberg A., 

2000, ApJ, 542, 235 
Laing R.A., Riley J.M., Longair M.S., 1983, MNRAS, 204, 

151 

Lampton M., Margon B., Bowyer S., 1976, ApJ, 208, 177 
Leahy J.P., Williams A.G., 1984, MNRAS, 210, 929 
Liu R., Pooley G., Riley J.M., 1992, MNRAS, 257, 545 
Longair M.S., Seldner M., 1979, MNRAS, 189, 433 
Machalski J., Chyzy K.T., Jamrozy M., 2004a, AcA, 54, 
249 

Machalski J., Chyzy K.T., Jamrozy M., 2004b, AcA, 54, 
391 

Machalski J., Chyzy K.T., Stawarz L., Koziel. D., 2007, 

A&A, 462, 43 
Manolakou K., Kirk J.G., [MK], 2002, A&A, 391, 127 
Markoff S., Falcke H., Fender R., 2001, A&A, 372, 25 
Martmez-Sansigre A., Rawlings S., 2011, MNRAS, 414, 

1937 

Massaro F., Ajello M., 2011, ApJ, 729, 12 

McGilchrist M.M, Baldwin J.E., Riley J.M., Titterington 

D.J., Waldram E.M., Warner P.J., 1990, MNRAS, 246, 

110 

McLure R.J., Dunlop J.S., 2004, MNRAS, 352, 1390 

McNamara B.R., Wise M., Nulsen P.E.J. , David L.P., 
Sarazin C.L., Bautz M., Markevitch M., Vikhlinin A., For- 
man W.R., Jones C, Harris D.E., 2000, ApJ, 534, 135 

Meier D.L., Edgington S., Godon P., Payne D.G.,Lind 
K.R., 1997, Nature, 388, 350 

Meier D.L., 1999, ApJ, 522, 753 

Meier D.L., 2002, NewAR, 46, 247 

Meisenheimer K., Roser H.-J., Hiltner P.R., Yates M.G., 
Longair M.S., Chini R., Perley R.A., 1989, A&A, 219, 63 
Meli A., Mastichiadis A., 2008, NIMPA, 588, 193 
Merloni A., Heinz S., 2008, MNRAS, 388, 1011 
Mullin, L.M., Riley J.M., Hardcastle M.J., 2008, MNRAS, 
390, 595 

Natarajan P., Treister E., 2009, MNRAS, 393, 838 
O'Dea CP, Daly R.A., Kharb P., Freeman K.A., Baum 

S.A., 2009, ApJ, 494, 471 
Penrose R., Floyd G.R., 1971, Nature Physical Science, 

229, 177 

Poggianti B.M., De Lucia G., Varela J., Aragon-Salamanca 
A., Finn R., Desai V., von der Linden A., White S.D.M., 
2010, MNRAS, 405, 995 



Press W.H., Teukolsky S.A., Vetterling W.T., Flannery 
B.P., 1992, 'Numerical recipes in C. The art of scientific 
computing', Cambridge: University Press, 2nd ed. 
Rawlings S., Jarvis M.J.., 2004, MNRAS, 355, 9 
Rawlings S., Saunders R., 1991, Nature, 349, 138 
Rees M.J., Begelman M.C., Blandford R.D., Phinney E.S., 

1982, Nature, 295, 17 
Rosati P., Tozzi P., Ettori S., Mainieri V., Demarco R., 
Stanford S.A., Lidman C, Nonino M., Borgani S., Delia 
Ceca R., Eisenhardt P., Holden B.P., Norman C, 2004, 
AJ, 127, 230 

Sanderson A.J.R., O'Sullivan E., Ponman T.J., 2009, MN- 
RAS, 395, 764 

Saripalli L., Subrahmanyan R., Hunstead R.W., 1994, MN- 
RAS, 269, 37 

Scannapieco E., Silk J., Bouwens R., 2005, ApJ, 635, 13 

Schechter P., 1976, ApJ, 203, 297 

Scheuer P.A.G., 1974, MNRAS, 166, 513 

Scheuer P.A.G., 1995, MNRAS, 277, 331 

Scheuer P.A.G., 1996, ASPC, 100, 1 

Schmidt M., 1968, ApJ, 151, 393 

Schoenmakers A. P., Mack K.-H., Lara L., Rottgering 
H.J. A., de Bruyn A.G., van der Laan H., Giovanni G., 
1998, A&A, 336, 455 
Shabala S.S., Kaviraj S., Silk J., 2011, MNRAS, 413, 2815 
Shankar F., Weinberg D.H., Miralda-Escude J., 2009, ApJ, 
690, 20 

Shklovskii I.S., 1963, SvA, 6, 465 

Sikora M., Begelman M.C., Madejski G., Lasota J.-P.. 2005, 

ApJ, 625, 72 
Silk J., 2005, MNRAS, 364, 1337 
Singal A.K., 1993, MNRAS, 263, 139 
Spitkovsky A., 2008, ApJ, 682, 5 
Steinhardt C.L., Elvis M., 2010, MNRAS, 402, 2637 
Steinhardt C.L., Elvis M., 2011, MNRAS, 410, 201 
Ueda Y., Akiyama M., Ohta K., Miyaji T., 2003, ApJ, 598, 

886 

van Ojik R., Rottgering H.J.A., Carilli, C.L., Miley G.K., 
Bremer M.N., Macchetto F., 1996, A&A, 313, 25 

Wang Y., Kaiser C.R., 2008, MNRAS, 388, 677 

Wang Y., Knigge C, Croston J.H., Pavlovski G., 2011, MN- 
RAS, in press, [arXivll08.0558W] 

Wardle J.F.C., Homan D.C., Ojha R., Roberts D.H., 1998, 
Nature, 395, 457 

Willott C.J., Rawlings S., Blundell K.M., Lacy M., 1999, 
MNRAS, 309, 1017 

Willott C.J., Rawlings S., Blundell K.M., Lacy M., Eales 
S.A., 2001, MNRAS, 322, 536 

Wold M., Lacy M., Lilje P.B., Serjeant S., 2000, MNRAS, 
316, 267 

Wolfe A.M., Burbidge G.R., 1970, ApJ, 161, 419 

Wu X.-B., Liu F.K., Zhang T.Z., 2002, A&A, 389, 742 

Zirbel E.L., 1997, ApJ, 476, 489 

This paper has been typeset from a T^X/ DT^X file prepared 
by the author. 



Fundamental properties of FR II radio galaxies via Monte Carlo simulations 31 

APPENDIX A: 

SUPPLEMENTARY 
MATERIAL 



32 A.D. Kapniska, P. Uttley & C.R. Kaiser 





Figure Al. 90 per cent joint confidence intervals for radio galaxy populations created with different assumptions on the standard 
deviation of the density log-normal distribution, where cr log 10 (p ) = 0.0 (i.e. p is a delta function, drawn in red), Ciogi (po) = 0.15 
(black), o"i oglo ( po ) = 0.50 (violet), ci O g 10 (p ) = 0.75 (green), and 0"i O g lo (p o ) = 1.0 (blue). Results shown for the redshift range < Z\ < 0.3 
and for Model S as described in §5.4.1 of the main paper. Findings are based on the 3CRR and BRL data fits. The case where the 
standard deviation is equal to ciog 10 (po) = 0-15 is the initial, standard case analysed in depth in §5 of the main article. 




Figure A2. 90 per cent joint confidence intervals for radio galaxy populations created with different assumptions on the standard 
deviation of the density log-normal distribution, where fi O g 10 (p ) = 0.0 (i.e. p is a delta function, drawn in red), "logio(po) = 0-1S (black), 
CT logio(po) = (violet), Ciog 10 (po) = (g reen )i an d °k>gio(po) = ^ (blue). Results shown for the redshift range 0.8 < 23 < 2.0 and 
for Model S as described in §5.4.1 of the main paper. Findings are based on the 3CRR and BRL data fits. The case where the standard 
deviation is equal to CT iog 10 (p ) = 0.15 is the initial, standard case analysed in depth in §5 of the main article. 
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Figure A3. 90 per cent joint confidence intervals for radio galaxy populations created with different assumptions on the standard 
deviation of the log-normal maximum lifetimes distribution of the sources, where fiog 10 (t max ) = 0-05 (black), "logm (Wax) = 0*3 (blue), 
and °"i O g 10 (t max ) = 0-6 (red). Results shown for the redshift range < 23 < 0.3 and for Model S as described in §5.4.2 of the main paper. 
Findings are based on the 3CRR and BRL data fits. The case where the standard deviation is equal to CT iog 10 (t max ) = 0-05 i s t ne initial, 
standard case analysed in depth in §5 of the main article. 
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Figure A4. 90 per cent joint confidence intervals for radio galaxy populations created with different assumptions on the standard 
deviation of the log-normal maximum lifetimes distribution of the sources, where fiog 10 (t max ) = 0-05 (black), ci O g 10 (t max ) = 0-3 (blue), 
and CT iog 10 (t max ) = 0-6 (red). Results shown for the redshift range 08 < 23 < 2.0 and for Model S as described in §5.4.2 of the main 
paper. Findings are based on the 3CRR and BRL data fits. The case where the standard deviation is equal to ""iog 10 (t max ) = 0-05 i s the 
initial, standard case analysed in depth in §5 of the main article. 
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Figure A5. 90 per cent joint confidence intervals for radio galaxy populations created with different assumptions of the maximum head 
advance speeds of jets: 0.8c (orange), 0.4c (black), 0.15c (blue), and 0.05c (violet) in units of the light speed c. Results shown for the 
redshift range < z\ < 0.3 and for Model S as described in §5.4.3of the main paper. Findings are based on the 3CRR and BRL data 
fits. The case when the maximum head advance speed is assumed to be 0.4c is the initial, standard case analysed in depth in §5 of the 
main article. 
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Figure A6. 90 per cent joint confidence intervals for radio galaxy populations created with different assumptions of the maximum head 
advance speeds of jets: 0.8c (orange), 0.4c (black), 0.15c (blue), and 0.05c (violet) in units of the light speed c. Results shown for the 
redshift range 0.8 < 23 < 2.0 and for Model S as described in §5.4.3 of the main paper. Findings are based on the 3CRR and BRL data 
fits. The case when the maximum head advance speed is assumed to be 0.4c is the initial, standard case analysed in depth in §5 of the 
main article. 
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Figure A7. 90 per cent joint confidence intervals for radio galaxy populations created with different assumptions on the particle injection 
index, where (i) m 6 [2; 3] (uniform distribution, drawn in black), (ii) m = 2.3 (single value, violet), and (iii) m = 2.4 and cr mo = 0.3 
(green). Results shown for the redshift range < z\ < 0.3 and for Model S as described in §5.4.4 of the main paper. Findings are based 
on the 3CRR and BRL data fits. The case where the particle injection index is assumed to be a uniform distribution, the case (i) above, 
is the initial, standard case analysed in depth in §5 of the main article. 
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Figure A8. 90 per cent joint confidence intervals for radio galaxy populations created with different assumptions on the particle injection 
index, where (i) m £ [2; 3] (uniform distribution, drawn in black), (ii) m = 2.3 (single value, violet), and (iii) m = 2.4 and cr mo = 0.3 
(green). Results shown for the redshift range 0.8 < 2:3 < 2.0 and for Model S as described in §5.4.4 of the main paper. Findings are based 
on the 3CRR and BRL data fits. The case where the particle injection index is assumed to be a uniform distribution, the case (i) above, 
is the initial, standard case analysed in depth in §5 of the main article. 
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Figure A9. 90 per cent joint confidence intervals for radio galaxy populations created with different assumptions on the standard 
deviation of the density log-normal distribution, where (i) 7 m i n = 1 and 7 max = 10 10 (drawn in black), (ii) 7 m ; n = 10 2 and 7 max = 10 10 
(magenta), (iii) 7 m i n = 10 4 and 7 max = 10 10 (violet), and(iv) 7 m i n = 1 and 7 ma x = 10 5 (green). Results shown for the redshift range 
< Zi < 0.3 and for Model S as described in §5.4.4 of the main paper. Findings are based on the 3CRR and BRL data fits. The case 
where the Lorentz factors are set to 7 m i n = 1 and 7 max = 10 10 is the initial, standard case analysed in depth in §5 of the main article. 





Figure A10. 90 per cent joint confidence intervals for radio galaxy populations created with different assumptions on the standard 
deviation of the density log-normal distribution, where (i) 7 m i n = 1 and 7 max = 10 10 (drawn in black), (ii) 7 m i n = 10 2 and 7 max = 10 10 
(magenta), (iii) 7 m i n = 10 4 and 7 max = 10 10 (violet), and(iv) 7 m i n = 1 and 7 max = 10 5 (green). Results shown for the redshift range 
0.8 < 23 < 2.0 and for Model S as described in §5.4.4 of the main paper. Findings are based on the 3CRR and BRL data fits. The case 
where the Lorentz factors are set to 7 m i n = 1 and 7 max = 10 lf) is the initial, standard case analysed in depth in §5 of the main article. 
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Figure All. 90 per cent joint confidence intervals for radio galaxy populations created with different assumptions on the particle content 
of the jet, where k' = (lightweight jets, black), k' £ [0; 10] (modest inclusion of protons in the jet, light green), k' € [0; 100] (heavy jets, 
blue), and k' = 100 (heavy jets, orange). Results shown for the redshift range < z\ < 0.3 and for Model S as described in §5.4.5 of the 
main paper. Findings are based on the 3CRR and BRL data fits. The case where k' = is the initial, standard case analysed in depth 
in §5 of the main article. 




Figure A12. 90 per cent joint confidence intervals for radio galaxy populations created with different assumptions on the particle content 
of the jet, where k 1 = (lightweight jets, black), k' £ [0; 10] (modest inclusion of protons in the jet, light green), k 1 S [0; 100] (heavy 
jets, blue), and k' = 100 (heavy jets, orange). Results shown for the redshift range 0.8 < z\ < 2.0 and for Model S as described in §5.4.5 
of the main paper. Findings are based on the 3CRR and BRL data fits. The case where k' = is the initial, standard case analysed in 
depth in §5 of the main article. 
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Figure A13. 90 per cent joint confidence intervals for radio galaxy populations created with various assumptions on the minimum 
energy of the initial particle distribution and the particle content of the jet. The following cases are plotted: (i) k' = 0, 7 m i n = 1 and 
7max = 10 10 (solid, black, the default case), (ii) k' £ (0, 100), y min = 10 2 and 7 max = 10 10 (solid, red), (iii) k' = 100, y min = 10 2 and 
7max = 10 10 (solid, orange), (iv) k' G (0, 100), 7min = 10 4 and 7max = 10 10 (solid, green), and (v) k' = 100, 7min = 10 4 and 7max = 10 10 
(solid, blue). For reference plotted are the cases where (a) k r = 100, 7 m i n = 1 and 7 max = 10 10 (dotted, dark grey, see also Figure lAHJl . 
and (b) k' = 0, 7 m in = 10 4 and 7 ma x = 10 10 (dotted, light grey, see also Figure \K9\i . Results shown for the redshift range < z\ < 0.3 
and for Model S as described in §5.4.5 of the main paper. Findings are based on the 3CRR and BRL data fits. 
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Figure A14. 90 per cent joint confidence intervals for radio galaxy populations created with various assumptions on the minimum 
energy of the initial particle distribution and the particle content of the jet. The following cases are plotted: (i) k' = 0, 7 m i n = 1 and 
7max = 10 10 (solid, black, the default case), (ii) k' £ (0, 100), y min = 10 2 and 7 max = 10 10 (solid, red), (iii) k' = 100, y min = 10 2 and 
7max = 10 10 (solid, orange), (iv) k' G (0, 100), 7min = 10 4 and 7max = 10 10 (solid, green), and (v) k' = 100, 7min = 10 4 and 7max = 10 10 
(solid, blue). For reference plotted are the cases where (a) k r = 100, 7 m i n = 1 and 7 max = 10 10 (dotted, dark grey, see also Figure |A12|| . 
and (b) k' = 0, 7 m in = 10 4 and 7 max = 10 10 (dotted, light grey, see also Figure lAlOI I, Results shown for the redshift range 0.8 < 23 < 2.0 
and for Model S as described in §5.4.5 of the main paper. Findings are based on the 3CRR and BRL data fits. 



